Magnetomicrometric Advances in Robotic Control

ABSTRACT

Systems and methods relating to magnetomicrometry and magnetic-target-based mechanomyography are provided. A method of detecting muscle activation, includes with a magnetic field sensor, detecting lateral vibration of a target implanted at a muscle or a tendon and estimating a level of muscle activation based on the detected lateral vibration. The target comprises a magnetic material.

RELATED APPLICATION(S)

This application claims the benefit of U.S. Provisional Application No.63/104,942, filed on Oct. 23, 2020. The entire teachings of the aboveapplication are incorporated herein by reference.

BACKGROUND

Target tracking is useful across a wide range of disciplines and atvarious scales, including, for example, tissue monitoring, beam bending,fluid dynamics, human-computer interaction, and traffic.

There are many methods to track visually-obscured objects.Magnetic-target tracking has the advantages of being low-cost, portable,and safe. However, current magnet tracking technologies are slow,precluding high-speed real-time magnetic-target tracking. There exists aneed for improved magnet tracking methods, as can be used, for example,in tissue monitoring applications.

SUMMARY

Methods and devices are disclosed that relate generally to magnetictargets and associated sensors and processes for use, for example, inmagnetomicrometry and mechanomyography.

Implantable Targets and Insertion Devices

An implantable target for magnet tracking includes a base structure anda shell structure. The base structure comprises a magnetic ormagnetizable material, and the shell structure comprises layers ofnickel, copper, gold, and parylene C.

An insertion device for an implantable target includes a cannula, acartridge, and a pushrod. The cartridge is configured to position animplantable target at the cannula, and the pushrod is receivable in thecannula and configured to push the implantable target from the cartridgethrough the cannula for delivery to a tissue. The cannula and thepushrod comprise a nonmagnetic material.

An insertion system can include a plurality of implantable targets andan insertion device. The insertion system can provide for delivery ofthe implantable targets.

Muscle Activation Detection

A method of detecting muscle activation includes, with a magnetic fieldsensor, detecting lateral vibration of a magnetic target implanted at amuscle or tendon and estimating a level of muscle activation based onthe detected vibration.

A system for detecting muscle activation includes a magnetic fieldsensor configured to detect lateral vibration of at least one targetimplanted at a muscle or a tendon, the at least one target comprising amagnetic material. The system further includes a controller configuredto estimate a level of muscle activation based on the detected lateralvibration.

Physiological Parameter Estimation

A method of estimating a physiological parameter of a muscle or tendon,or a muscle-tendon unit, includes, with a magnetic field sensor,detecting vibration of a magnetic target implanted at a muscle or tendonand estimating at least one of a muscle force and tendon force based onthe detected vibration.

A system for estimating a physiological parameter of a muscle or tendon,or a muscle-tendon unit, includes a magnetic field sensor configured todetect vibration of at least one target implanted at a muscle or tendon,the at least one target comprising a magnetic material. The systemfurther includes a controller configured to estimate at least one of amuscle force and tendon force based on the detected vibration.

Portable Goniometry

A method of monitoring biomechanical motion includes disposing at leastone target on a subject at a location associated with a joint of thesubject and, with a magnetic field sensor array, detecting a change instate of the at least one target relative to the magnetic field sensorarray, another target disposed on the subject, or a combination thereof.The method further includes determining a state of the joint based onthe detected change in state of the at least one target.

Improvements to Magnetic Object Tracking

A method for determining one or more of three sensor position parametersand three sensor orientation parameters for each of the sensors in asensor array includes placing at least one target in at least one knownlocation relative to a sensor array, whereby a signal from the at leastone target at the sensors is detected, and recording at least onemeasurement of the signal at each of the sensors for each placement ofthe one or more targets. The method further includes estimating one ormore parameters from the group consisting of x-position, y-position,z-position, yaw, pitch, and roll, of each of the sensors and estimatingany unknown state parameters of the at least one target. A constantvalue for a magnetic dipole weight state parameter of the at least onetarget for each of the measurements is provided, and predicted values ofthe signal at each of the sensors for each of the measurements given theone or more estimated sensor parameters, the estimated target stateparameters, and the provided constant magnetic dipole weight stateparameter are calculated. A prediction error in the predicted values ofthe signal with reference to the values of the signals detected at thesensors is computed. A prediction error Jacobian matrix is calculated byanalytically computing elements of the prediction error Jacobian matrixwith respect to the estimated parameters of the sensors for eachmeasurement and, from the prediction error and the prediction errorJacobian matrix, a state of the parameters of the sensors is determined.

A method of tracking one or more permanent magnets includes detecting asignal from each of the one or more permanent magnets, calculating ananalytically-derived Hessian matrix with respect to the detectedsignals, and determining a state of each of the one or more permanentmagnets based on the calculated Hessian matrix.

A system comprises at least two sensor arrays, each sensor arrayconfigured to detect a state of at least one magnetic target at atissue, and at least one position sensor associated with at least one ofthe at least two sensor arrays and configured to detect a position andorientation of the associated sensor array relative to the other of theat least two sensor arrays.

Calibration Improvements for Magnetic Field Sensor Arrays

A method of calibrating a magnetic field sensor array that comprises aplurality of magnets includes three-dimensionally rotating a magneticfield sensor array in a uniform magnetic field and recording data fromeach of the plurality of sensors. The method further includescalculating a non-rotating transformation of the recorded data using anellipsoid fit of each of the plurality of sensors and calculating ascaling factor of each of the plurality of sensors based on relativedimensions of the ellipsoid fit. A transformed, scaled dataset for eachof the plurality of sensors is obtained through application of thecalculated non-rotating transformation and calculated scaling factor.The obtained transformed, scaled datasets are rotated to align, therebyproviding for a determination of a relative orientation of each of theplurality of sensors to calibrate the magnetic field sensor array.

Wearable Shielding and Non-Planar Sensing Arrays

A wearable shielding assembly comprises an array of sensors configuredto detect a state change of at least one magnetic target implanted at atissue, a wearable receptacle within which the array of sensors isdisposed, and a geometrically-reconfigurable material disposed about orintegral with the wearable receptacle. The geometrically-reconfigurablematerials provide magnetic shielding to the array of sensors.

A method of assembling a non-planar sensing array, such as for aprosthetic device, includes fabricating a plurality of sensors on aflexible circuit board and affixing the flexible circuit board is to arigid substrate.

BRIEF DESCRIPTION OF THE DRAWINGS

The patent or application file contains at least one drawing executed incolor. Copies of this patent or patent application publication withcolor drawing(s) will be provided by the Office upon request and paymentof the necessary fee.

The foregoing will be apparent from the following more particulardescription of example embodiments, as illustrated in the accompanyingdrawings in which like reference characters refer to the same partsthroughout the different views. The drawings are not necessarily toscale, emphasis instead being placed upon illustrating embodiments.

FIG. 1 is a schematic of an example magnetomicrometry andmechanomyography system that includes magnetic targets that passivelytransmit position information through the tissue to magnetic fieldsensors located external to the body.

FIG. 2 is a schematic of a wearable robot controlled, at least in part,by a magnetomicrometry system.

FIGS. 3A and 3B are photos of an example magnetic field sensing arrayand which was used in experimentation. Two 6×8 magnetic field sensorgrids were custom designed and held together using a 3D-printed fixtureand nylon nuts and bolts (FIG. 3A). The sensor positioning fixture wasalso used to house a custom adapter board (FIG. 3B) which connected themicrocontroller to the sensors via flat flexible cables.

FIG. 4 is a schematic of an experimental set up for comparison ofprototype magnetomicrometry devices with fluoromicrometry. A motor wasused to apply a mechanical frequency sweep from 0.7 to 7 Hz, with aspring to provide an opposing force. A laptop computer and a magneticfield sensor array mounted external to the turkey's leg were used totrack the distance between the magnetic beads (magnetomicrometry), andfluoromicrometry (X-ray stereo videofluoroscopy) was used to obtain acomparison measurement of the distance between the beads.

FIGS. 5A-D are graphs of measured distance over time for first (FIG.5A), second (FIG. 5B), third (FIG. 5C), and fourth (FIG. 5D) trialsperformed with a turkey (Turkey A). In each graph, magnetomicrometry(blue, line A) against fluoromicrometry (orange, line B) is plotted,with absolute difference (green, line C), for the frequency sweepperformed on the gastrocnemius of Turkey A. All plots correspond totrials from the same leg (right side).

FIGS. 6A-6F are graphs of measured distance over time for first (FIG.6A), second (FIG. 6B), and third (FIG. 6C) trials performed on the leftleg of another turkey (Turkey B), and first (FIG. 6D), second (FIG. 6E),and third (FIG. 6F) trials performed on the right leg of the turkey. Ineach graph, magnetomicrometry (blue, line A) against fluoromicrometry(orange, line B) is plotted, with absolute difference (green, line C),for the frequency sweep performed on the gastrocnemius of Turkey B.

FIGS. 7A-7F are graphs of measured distance over time for first (FIG.7A), second (FIG. 7B), and third (FIG. 7C) trials performed on the leftleg of yet another turkey (Turkey C), and first (FIG. 7D), second (FIG.7E), and third (FIG. 7F) trials performed on the right leg of theturkey. In each graph, magnetomicrometry (blue, line A) againstfluoromicrometry (orange, line B) is plotted, with absolute difference(green, line C), for the frequency sweep performed on the gastrocnemiusof Turkey C.

FIGS. 8A-8F are graphs of measured distance over time for first (FIG.8A), second (FIG. 8B), and third (FIG. 8C) trials performed on the leftleg of yet another turkey (Turkey D), and first (FIG. 8D), second (FIG.8E), and third (FIG. 8F) trials performed on the right leg of theturkey. In each graph, magnetomicrometry (blue, line A) againstfluoromicrometry (orange, line B) is plotted, with absolute difference(green, line C), for the frequency sweep performed on the gastrocnemiusof Turkey D.

FIG. 9 is a histogram of differences between the magnetomicrometry andfluoromicrometry gastrocnemius frequency sweep measurements, inmicrometers. Different histogram colors correspond to all four trialswith Turkey A right leg (green), all three trials with Turkey B left legand right legs (purple and red, respectively), all three trials withTurkey C left and right legs (blue and orange, respectively), and allthree trials with Turkey D left and right legs (cyan and brown,respectively). The mean difference and standard deviation of thedifference, in micrometers, is shown to the right of the histogram.

FIG. 10 is a grid of histograms of magnetomicrometry andfluoromicrometry static measurements. Fluoromicrometry (orange) andmagnetomicrometry (blue) measurements were taken while magnets wereplaced at separation distances of approximately 24, 40, 56, and 72 mm ina LEGO block (vertical dashed gray lines in each column), with thesensing array at various heights above the magnets (sensor proximity,shown to the right of the plots.

FIG. 11 is a grid of histograms of magnetomicrometry andfluoromicrometry, static measurements. Fluoromicrometry (orange) andmagnetomicrometry (blue) measurements were taken while magnets wereplaced at separation distances of approximately 24, 40, 56, and 72 mm ina LEGO block (horizontal dashed gray lines in each row), with thesensing array at various heights above the magnets. In contrast withFIG. 10 , the vertical axis shows the separation distances, and thesensor proximity is shown at the top of the plots. Magnetomicrometrydata is plotted for one second before and one second after the start andend of the fluoromicrometry measurements, respectively. Note that someof the fluoromicrometry measurements also exhibit a warm-up time lastinga small fraction of a second.

FIG. 12 is a plot of histograms of single-board magnetomicrometry staticmeasurements. Magnetomicrometry measurements (blue) were taken whilemagnets were placed at a separation distances of approximately 24 mm ina LEGO block (indicated by the vertical dashed gray line), with thesensing array at various heights above the magnets (sensor proximity),shown to the right of the plots.

FIG. 13 is a graph of separation distance over time for long-termimplant stability of 3mm-diameter magnet pairs against migration inmuscle. Pairs of 3mm-diameter magnets were implanted with variousseparation distances into the gastrocnemius (blue lines, labelled “G”)and iliotibialis cranialis (brown lines, labelled “IC”) muscles of fourturkeys. Separation distances were monitored over time via computedtomography scans. Note that there is a cutoff point (around 20 mm forthe 3-mm-diameter magnets used) where magnets should not be implantedany closer to one another, to ensure stability against migration.

FIG. 14 is a graph of force between 3-mm-diameter magnets as a functionof magnet separation distance. Force between a pair of 3-mm-diametermagnets (shown in blue) at various separation distances, assumingalignment along the line of separation between the magnets. Blue shadingrepresents the forces calculating using the magnetic dipoles strengthsplus and minus their respective measurement standard deviations. Forreference, the force due to gravity on each magnet is indicated by anorange horizontal line.

FIG. 15 is a schematic of an example magnetic target.

FIG. 16 is a flow diagram of an example method of magnetic-target-basedmechanomyography.

FIG. 17 is a diagram illustrating an example control scheme usingmagnetic-target-based mechanomyography and magnetomicrometry for controlof a robotic prosthesis.

FIG. 18 is a schematic of a portable goniometer.

FIG. 19 is a schematic of a magnetic bead insertion device.

FIG. 20 is a schematic of a phantom magnet below a circular sensor arrayconfiguration as a visual of a proof performed. A circular sensor arrayis shown measuring the magnetic field from a north-side-up magnetcentered at a height of approximately 0.707 times the sensor circleradius. A north-side-down magnet below the array (outlined in dashedblack) creates the same magnetic fields at each of the sensors, and thusis an alternative solution, or “phantom” magnet.

FIG. 21 is a flow diagram illustrating an example method of performing arefined dipole measurement.

FIG. 22 is a schematic of a magnetic field sensor and rotation in aspatially-uniform ambient magnetic field for use in accounting for hardand soft iron effects.

FIGS. 23A-C are plots of calibration data after ellipsoid-to-spheretransformation. Each plot shows the three-dimensional calibration data,after ellipsoid-to-sphere transformation, from a different axis, inparticular, as viewed from the positive z-direction (FIG. 23A), positivex-direction (FIG. 23B), and negative y-direction (FIG. 23C). Distinctcolors/shading of the dots correspond to the different sensors in thearray. Note how the transformation results in data streams which areapproximately aligned, and how the data is approximately spherical andcentered at the origin.

FIG. 24 is a flow diagram of a calibration process including relativesensor rotation.

FIG. 25 is a graph illustrating measured ambient magnetic fielddemonstrating the magnetic field effects of a passing train. Labels onthe right correspond to the x, y, and z axes from top to bottom. Datawas sampled at 300 Hz, so this measurement covers about 40 seconds ofmagnetic field measurement, with the train passing event lasting about10-15 seconds.

FIG. 26 is a schematic of an example wearable shielding assemblyincluding geometrically-reconfigurable material.

FIG. 27 is a schematic of another example wearable shielding assemblyincluding geometrically-reconfigurable material.

FIG. 28 is a flow diagram representing common subexpression eliminationfor Jacobian matrix elements. The color represents the number offloating point operations required to get to each variable, in theperceptionally uniform viridis colormap, incrementing from bright yellowto dark purple. To appreciate the efficiency gained from commonsubexpression elimination, note, for instance, how x ⁻²ij, y ⁻²ij, z⁻²ij are each only computed once, as opposed to the 212 total times thatthey would be computed using Eqns. 11a-r in the Taylor Reference.

FIG. 29 is a flow diagram representing common subexpression eliminationfor orientation parameters. The color represents the number of floatingpoint operations required to get to each variable, in the perceptuallyuniform viridis colormap, incrementing from bright yellow to dark purple(though it is noted that sines and cosine generally requiresignificantly more time to implement than, for instance, addition ormultiplication). All of the output variables from this flow diagram canbe used in the computation of the Jacobian matrix elements, whereas onlyc₀, c₁ and c₂ are used in the calculation of the magnetic fieldprediction.

FIG. 30 is a flow diagram representing common subexpression eliminationfor magnetic field predication calculation. The color represents thenumber of floating point operations required to get to each variable, inthe perceptionally uniform viridis colormap, incrementing from brightyellow to dark purple. Note that c₀, c₁ and c₂ are calculated as shownin FIG. 29 .

FIG. 31 is a schematic of a system comprising multiple sensing arraysfor which positions, or poses, can be tracked.

FIG. 32 is a flow diagram of an example method of monitoring poses of amultiple sensing array system.

FIG. 33 is a flow diagram of an example method of ranking sensors of amulti-sensor system for a target.

FIG. 34 is a schematic of a non-planar sensing array.

FIG. 35 is a Jacobian matrix constructed with a single magnet measuredby N sensors over K timesteps.

FIG. 36 is a sensor position calibration matrix.

DETAILED DESCRIPTION

A description of example embodiments follows.

Magnetomicrometry is a technology that tracks visually-obscured magneticbeads implanted within or on biological tissue. For example,magnetomicrometry methods and devices described herein can be applied tomonitor in-vivo tissue length and speed within freely moving animals andhumans.

The methods and devices described herein can be implemented in varioustechnologies, such as for the treatment of limb pathology resulting fromdisease or traumatic injury and for human augmentation to enhance humanphysicality beyond normal physiological limits. In the realm ofpermanent assistance devices, the methods and devices described hereincan be implemented in technologies for the preservation ofpost-amputation function in the residuum for the case of limbamputation, or for the restoration of natural muscle control function inparalyzed or weakened limbs due to age-related degeneration, spinal cordinjury, or other neuromuscular pathologies.

Descriptions of systems and methods relating to tracking of magneticobjects can be found in WO2019/074950, “Method for Neuromechanical andNeuroelectromagnetic Mitigation of Limb Pathology;” the entire teachingsof which are incorporated herein by reference.

Additional descriptions of systems and methods relating to tracking ofmultiple targets, such as permanent magnets, can be found in thefollowing publication: Cameron R Taylor, Haley G Abramson, and Hugh MHerr. Low-latency tracking of multiple permanent magnets. IEEE SensorsJournal, 19(23):11458-11468, 2019 (the “Taylor Reference”); the entireteachings of which are incorporated herein by reference. In someaspects, the features described herein present improvements orextensions of the methods and devices described in WO2019/074950 and inthe Taylor Reference.

It should be noted that any of the features described herein can becombined with one another and with the features described in the notedpublications, in any combination, in a method or system. For example,the monitoring of poses of a multiple sensing array system, as shown inFIG. 32 , can be combined with mechanomyographic methods, as shown inFIG. 16 . While some examples herein refer to single amputation levelsin single extremities, the methods and devices described can be appliedto other limbs and amputation levels.

The terms “magnetic bead” and “magnetic target” are used interchangeablyherein. Where the term “magnetic bead” is used, it should be understoodthat the magnetic target need not be of a spherical geometry, andtargets of other shapes are possible and included within the meaning ofthe term “magnetic bead,” including, for example, cube magnets

As defined herein, a “state” of a target, or targets, includes at leastone of any of the members of the group consisting of the position,orientation, and strength of the target or targets. Specifically, the“state of the targets relative to each other” includes at least one ofany of the members of the group consisting of the relative positions ofthe targets, the distance between the targets, and the relativeorientation between the targets.

There are many methods to track visually-obscured objects.Magnetic-target tracking has the advantages of being low-cost, portable,and safe. However, current magnet tracking technologies are slow,precluding high-speed real-time magnetic-target tracking. This is due tothe mathematics of magnet tracking, whereby magnet positions aretraditionally determined via numerical optimization, which can sufferfrom instability and significant delays. Improved methods and systemsfor tracking one or more magnets with high speed and accuracy areprovided. Additionally, validation of such methods is provided throughdemonstrations of real-time muscle length tracking.

The methods and systems described herein can provide for high-speed,real-time, multiple-magnetic-target tracking. Such methods and systemscan use an analytic gradient of a magnetic field prediction error.Magnetic disturbances can be compensated for in real time using asimpler, more portable strategy than currently-published disturbancecompensation methods. Validating the method in a physical system againststate-of-the-art motion capture, increased maximum bandwidths of 336%,525%, 635%, and 773% were demonstrated for the simultaneous tracking of1, 2, 3, and 4 magnets, respectively, with tracking accuracy comparableto state-of-the-art magnet tracking, as further described herein.

Using pairs of implanted magnetic beads to wirelessly track musclelength and speed, a mechanical frequency sweep was applied to an in-vivoturkey gastrocnemius muscle and submillimeter agreement betweenmagnetic-bead-derived real-time muscle length measurements and stereoX-ray videofluoroscopy was found. Longitudinal data using computedtomography was collected, and a minimum magnetic-bead separationdistance in muscle of approximately two centimeters was found.

Magnetomicrometry can provide for peripheral nervous system control ofwearable robots via real-time tracking of muscle lengths and speeds, aswell as for the in-vivo tracking of biological tissues to elucidatebiomechanical principles of animal and human movement.

Magnetomicrometry can be used to deliver intuitive, skillful controlover bionic prostheses. Magnetomicrometry was developed to meet animmense need in prosthetic control. In the United States alone, thereare likely 2 to 3 million persons living with the loss of a limb [25].However, the commercially-available methods for controlling prostheticlimbs lag behind current robotic prosthesis technologies. On one side ofthe spectrum are fully non-invasive technologies, such aselectromyography, ultrasound, and mechanomyography, all of which resideoutside the body, but can have poor, unstable signal quality. On theother side of the spectrum are highly-invasive technologies, such assonomicrometry, micro-channel arrays, and intrafascicular electrodes.These highly invasive technologies provide improved signal quality, butare expensive to implement and require delicate surgery and carefulpost-operative care. The field is currently missing a low-footprint,minimally-invasive tissue interface that can accurately monitor muscleactions, as well as muscle activation.

Systems and methods relating to magnetomicrometry, a new technologydeveloped for tracking tissue lengths, speeds and properties, arepresented. Magnetomicrometry generally involves implanted magneticbead(s), which can be used to wirelessly track tissue lengths and tissuestates, making it the first ever minimally-invasive real-time muscletracking technology. Magnetomicrometry systems and methods can also beimplemented with simple tracking hardware, making it economical,compact, and portable. With example prototype systems, real-time musclelength tracking with sub-millimeter accuracy and with precision towithin a tenth of a millimeter has been demonstrated, as describedfurther herein.

Significant progress has been made recently in the field of prostheticcontrol and in developing robotic prostheses. State-of-the-art devices,such as, for example, the robotic bebionic® hand (Ottobock, Austin, TX),have many degrees of freedom and are strong, durable tools. However,these devices fall short in their ability to provide full, intuitive,human-like control to the amputee.

The commercially-available method for controlling these advancedprostheses uses electrodes—bare metal conductors placed on the surfaceof the residual limb. These electrodes measure a signal called muscleactivation, and this activation signal is used to control the speed ofthe prosthesis. Any change in electrode placement results in requiredrecalibration, and even moisture changes throughout the day can degradethe signal. Though our lab, the Biomechatronics Group, has demonstratedsignificant improvements to this method, called electromyography (EMG),such an approach to prosthesis control has two fundamental limitations.

To accurately map muscle activation to muscle force, both the musclelength and the muscle velocity should be known [8]. A problem arises inthat there currently does not exist a technology for measuring theseparameters in humans in real time, and it is not possible to accuratelydetermine muscle force, length, or speed from EMG information. Thus,using EMG alone, delivering complete biomimetic control over theprosthesis is difficult. With sensing muscle length and velocity,accurately determining muscle force commands from muscle activation canbe provided, or, further, force control can be skipped and directposition control of the prosthesis can be provided. An objective of theprovided methods and systems is to provide for full, intuitive,human-like control of a prosthesis by improving the fidelity of forceand position control via muscle length and velocity measurements, whileeliminating or reducing perceptible control delays.

Magnetomicrometry enables wireless sensing of tissue lengths. Forexample, a magnetomicrometry approach can use pairs of magnetic beadsimplanted in tissue, which passively transmit position information tomagnetic field sensors external to the body. This technology wirelesslysenses muscle lengths and speeds in a low-cost, compact, portable,passive, and safe manner.

Examples of low-latency tracking of multiple permanent magnetics withdisturbance compensation are provided in the Taylor Reference [17], theentire teachings of which are incorporated herein.

Section 1 herein describes the results of preclinical work in turkeyswith prototype methods and systems for low-latency magnet tracking.Features relating to biocompatibility, accuracy, and long-term implantstability are provided. Section 2 provides for further description ofthe challenges of magnet tracking, and Section 3 provides for methodsfor magnetic dipole strength measurement to provide for improved magnettracking.

Section 4 provides for improvements to calibration of magnetic trackingsystems and methods. Section 5 provides for how magnetic fielddisturbance techniques can be extended to more general disturbances.Section 6 provides for methods for common subexpression eliminationtowards time delay reduction of magnet tracking methods and trackingoptimization. Lastly, Section 7 provides for further generalimprovements to magnet tracking methods and systems.

1. Minimally-Invasive Muscle Tracking Using Permanent Magnets

Low-footprint, minimally invasive tissue interfaces are provided thatcan accurately monitor muscle actions and that can overcome the variouslimitations of currently-available technologies. For example,fluoromicrometry, which uses X-rays for high precision tissue lengthtracking, is wireless but is limited to short bursts due to the ionizingradiation used and requires an entire room and significant processingtime [1]. In contrast, sonomicrometry uses physically-tethered,implanted ultrasound crystals to yield high accuracy [14], but it isinvasive and not easily miniaturized. And, electromyography, whether viasurface electrodes or highly-invasive direct-to-nerve interfaces, onlysenses muscle activation, which, without muscle length and velocity,cannot be used alone to faithfully observe, understand, or utilizemuscle action [8].

Magnetomicrometry is a new technology for tracking tissue lengths andspeeds. Implanted magnetic beads are used to wirelessly track tissuelengths, making magnetomicrometry the first-ever minimally-invasivereal-time muscle tracking technology. Furthermore, magnetomicrometrysystems and methods can use simple tracking hardware, making suchsystems economical, compact, and portable. With this system, real-timemuscle length tracking with sub-millimeter accuracy and with precisionto within a tenth of a millimeter has been demonstrated. This newtechnology, magnetomicrometry, lays the groundwork for full, biomimeticcontrol of prostheses and exoskeletons and an increased understanding ofthe biomechanical principles of animal and human movement.

Researchers have previously proposed implanting single magnets intomuscles in an attempt to monitor movement [16], but thesingle-magnet-per-muscle approach is limited in various ways. Musclelength can be passively cycled by the motion of a joint, such as theelbow joint when engaged by an aggressive handshake, or actively cycledwhen flexed, such as when holding a glass of water. These two movementscan confound one another if trying to measure either of them via theaxial motion of a single magnet, and lateral magnet position can bechallenging to use with a prosthesis due to changing pressures from thesocket and surrounding tissues. Finally, any motion of the magneticfield sensors or residual limb relative to one another can translatedirectly into tracking error. These issues can be solved by the use ofmagnet pairs in each muscle.

In this new tracking strategy, pairs of magnetic beads are wirelesslytracked as described in the Taylor Reference, with lightweight magneticfield sensor arrays mounted to a prosthetic socket. The compact hardwareused makes this technology not only portable but alsominimally-invasive. Because magnetic fields are not affected bymaterials such as silicon, carbon fiber, or the human body, the magneticfield can pass undisturbed from the muscles to the sensors as if theseother materials were not present.

An example magnetomicrometry system 100 is shown in FIG. 1 . The system100 includes implantable targets 101, 102 (e.g., magnetic beads). Asillustrated, the targets 101, 102 are implanted into biological tissue,such as a muscle, e.g., a muscle 105, or a tendon. The targets 101, 102can passively transmit position information through tissue to magneticfield sensors 106 a-e disposed external to the body. The magnetic fieldsensors 106 a-e can be in operative arrangement with a controller 110 toprovide measurements for use with magnetomicrometry and/ormechanomyography methods, as described further herein.

An example of a magnetomicrometry system for prosthetic control is shownin FIG. 2 . The system 200 includes pairs of implantable targets 201,202, and, optionally, an unpaired implantable target 203 (as can be usedfor mechanomyography, as described later herein). Magnetic field sensingarrays 208 a-c are disposed external of the body. The magnetic fieldsensing arrays can each include one or more magnetic field sensors 206and, optionally, at least one sensor 212 configured to detect a positionand/or an orientation of the array (e.g., an accelerometer, an inertialmeasurement unit, etc.). Muscle length, velocity, and activation can bemeasured via detected movement of the implanted targets and used toimpart control over a prosthesis 210 (e.g., as illustrated, a robotichand). The magnetic field sensing arrays 208 a-c can be in operativearrangement with a controller 210 (which can be disposed within theprosthesis 210 or an interface 215) to provide measurements for use withmagnetomicrometry and/or mechanomyography methods, as described furtherherein.

An in-vivo turkey model was investigated to address biocompatibility,verify in-vivo tracking accuracy, and ensure long-term resistance of theimplants to migration. The study methods and obtained results arefurther described in Example 1 in the Exemplification section below.

FIGS. 3A and 3B illustrate an example magnetic field sensing array 250,as was used in the experimentation described in Example 1. The magneticfield sensing array 250 includes two magnetic field sensor grids 252,254, disposed within fixture 256. The sensor positioning fixture wasalso used to house a custom adapter board 258 which connected themicrocontroller 260 to the sensors via flat flexible cables.

1.1 MRI-Safe Sensing Components

As described in Example 1, MRI-safe capacitors were used in the sensorarrays of the prototype devices. Traditional magnet tracking compensatesfor “hard” and “soft” iron effects, where some components that arerigidly affixed relative the magnetic field sensors cause unwanteddisturbance by retaining magnetization or by warping the magnetic fieldseen by the magnetic field sensors. To combat this, the use of circuitcomponents which do not warp or retain the magnetic field as seen by themagnetic field sensors can be used in magnetomicrometry andmechnomyography systems and devices.

For example, MRI-safe sensing components can be used for electricalcomponents on the sensing array board. Specifically, the use of MRI-safecapacitors (e.g. those which contain little to no ferromagnetic orotherwise magnetizable material) can be used in magnetic field sensors.These MRI-safe capacitors can be on the same or opposite side (top orbottom) of the sensing board circuitry and may be shared between sensorsas needed.

This use of MRI-safe capacitors is especially useful when proximitybetween sensors is tight, because a tight sensor layout can cause thecapacitors to be positioned especially close to the sensors. If thecapacitors distort the field, this increased proximity enlarges thisdistortion. With MRI-safe capacitors, the field distortion can beeliminated or reduced such that is no longer a factor for considerationin the positioning of the capacitors and sensors.

1.2 Biocompatibility

An example of an implantable target for magnetic tracking and providingfor biocompatibility is shown in FIG. 15 . The target 300 includes abase structure 301 that comprises a magnetic or magnetizable material.For example, the base structure can be formed from a neodymium ironboron+dysprosium base material. The target 301 further includes a shellstructure 302 disposed about the base structure that includes layers ofnickel, copper, gold, and parylene C.

The shell structure can include layers arranged, from the basestructure, in order of nickel, copper, nickel, gold, and parylene C. Agold layer of the shell structure can have a thickness of at least about5 μm. A paralyene C layer of the shell structure can have a thickness ofat least about 25 μm.

An example shell structure is shown in FIG. 15 . In the illustratedexample, a neodymium iron boron+dysprosium base material is coated innickel, copper, and then 99.99% pure conventional nickel plating,followed by a layer of gold (e.g., following ASTM B488, Type III, CodeA, Class 5 (99.9% Pure Gold, ≤90 HK25 Hardness, at least 5 micrometersthick)), followed by a layer of AdPro adhesion promoter and a coating ofParylene C (e.g., at least 25 micrometers thick).

The gold layer of the shell structure can provide a backup in case ofParylene C failure, and, in addition, can increase the radio-opaquenessof the implant, enhancing the ability to see the implant via alternativetechnologies (e.g., static x-ray or CT scans). All of these coatingsteps can be performed with the magnetic beads unmagnetized, with amagnetization step at the very end of the manufacturing process, beforeor after sterilization.

1.3 Muscle Force Measurements: Implanted-Magnetic-Bead-BasedMechanomyography

As discussed above, magnetomicrometry is a tool for measuring musclelength and velocity in real time, and additional information can beneeded to discern muscle intents and actions. To do this,magnetomicrometry can be paired with classical strategies, such assurface or needle electromyography, intrafascicular electrodes, ormechanomyography, to determine muscle force. Alternatively, or inaddition, muscle force measurements can be obtained from implantedmagnetic targets.

Previous work recording vibrations from isolated muscles in saline bathshas revealed some underlying principles of mechanomyography, whereinmuscle activation can be monitored by acoustic vibrations at the skin'ssurface using a microphone [3]. These experiments showed that, whenactivated, each muscle exhibits a lateral vibration (sideways, like aguitar string) on the order of 20-150 Hz. Thus, because the entiremuscle vibrates laterally during activation at a frequency higher thanthat of the muscle flexion itself, the lateral vibration of the samemagnetic beads used for muscle length and velocity sensing can bemonitored. This lateral vibration frequency can then be used to estimatemuscle activation, and, when combined with the length and velocity ofthe muscle, muscle force can be estimated using a biophysical musclemodel (e.g., Hill Muscle Model, [2]). The detection and measurement oflateral vibrations can depend on the precision of magnetomicrometry andthe amplitude of the lateral vibrations. The magnetomicrometry systemdescribed in Example 1 appears to be capable of monitoring lateralvibrations as low as around 10 to 20 micrometers, which can besufficient for providing muscle force measurements.

A method of detecting muscle activation includes, with a magnetic fieldsensor (e.g., sensors 106 a-e, or sensor arrays 208 a-c), detectinglateral vibration of a target implanted at a muscle or tendon (e.g.,targets 101, 102, 201, 202, 203). The method further includes estimatinga level of muscle activation based on the detected lateral vibration.

For example, detecting lateral vibration can include detecting movementsof the target in a range of about 10 μm to about 20 μm. Alternatively orin addition, detecting lateral vibration can include detectingvibrational movement of the target at frequencies of greater than about10 Hz. The detection of lateral vibration can be of a single target ormultiple targets. For example, detecting lateral vibration can includedetecting lateral vibration of two or more targets disposed along anaxis (e.g., axis A, FIG. 1 ). Optionally, an accelerometer can beincluded to detect vibrational movement of the magnetic field sensorrelative to the target. In combination with magnetomicrometry todetermine a length and a velocity of the muscle, as described above, amuscle force can be estimated based on the estimated level of muscleactivation.

Mechanomyography, or acoustic myography, has been used previously forsensing of muscle activation via vibration sensors or acoustic sensors,but it has never been sensed previously through the tracking of one ormore magnetic beads. One or more magnetic beads implanted in muscle(such as magnetic beads used for magnetomicrometry) can be employed toconvey mechanomyographic signals as magnetic field information tomagnetic field sensors external to the muscle. Such information can beconveyed to a computer via a wired or wireless connection, the computerand sensors being powered by a battery or some other power source.Specifically, mechanomyographic signals can result from muscle vibrationalong a single lateral dimension during muscle flexion. The lateraldimension is generally perpendicular to the longitudinal axis of themuscle.

For example, one or more magnetic beads can be implanted in a muscle,and, during flexion, these one or more magnetic beads undergo vibrationalong a single lateral dimension. The intensity and frequency of the oneor more magnets' vibrations can be determined using the magnetic fieldinformation conveyed from the magnetic field sensors. Because thevibrations occur at a frequency higher than biomechanical movements, thelower frequencies (<10 Hz) can be filtered out and the higherfrequencies (>10 Hz) can be used to determine muscle activation. Thisactivation signal from mechanomyography can then be combined with musclelength and velocity signals from, for example, magnetomicrometry,providing a robust estimate of muscle force using a biophysical musclemodel (e.g., the Hill Muscle Model) capable of predicting muscle forcegiven inputs of muscle length, velocity and activation.

For example, each of one or more magnetic beads implanted in or ontissue can be monitored for position, and the three components of theposition of each magnet can be separately analyzed for frequency contentusing, for instance, a Fourier or wavelet transform. These frequencies,between different components and between different magnets, can then becombined by, for example, weighted averaging to arrive at a finalfrequency with which to calculate muscle activation.

As a further example, the positions of multiple beads can be averagedbefore analyzing the frequency content of the three position componentsof the signal. For example, the previous N samples (where N is aninteger) can be used to fit a simple linear regression in 3D-space usingthe three position components for each of the magnets or for theaveraged position, and the frequency content of the resultingone-dimensional spatial information can then be analyzed.

More than one magnetic bead can be used, and knowledge of the placementof the magnetic beads can be employed to improve sensing capability. Forexample, with one magnetic bead distal and one magnetic bead proximalwithin a muscle, a line between the distal and proximal magnetic beadpositions (e.g., a time-averaged line, or a regression line includingadditional beads) can be used to reduce the dimensions along which dataanalysis occurs by projecting the data from each bead onto a planeorthogonal to the line between the beads. The resulting two dimensionsof position data for each bead can then be analyzed for frequencycontent separately, or can be reduced to one dimension using, forexample, a simple linear regression through the previous N data points(where N is some integer), before being analyzed for frequency content.Alternatively, the position data from the multiple magnetic beads can beaveraged before or after projection into 2D-space. Regardless of themethod, a mean peak frequency can be output, and, optionally, anamplitude at this frequency. This single frequency, as well as itsamplitude, if calculated, can then provide for calculation of muscleactivation.

An example method 400 of magnetic-target-based mechanomyogrpahy is shownin FIG. 16 . The method uses, as an example, a pair of implantedmagnetic targets. The method includes, for each timestep, storing thethree-dimensional positions of the two, tracked magnetic beads (402). Avector is computed from a time-averaged position of one of the two beadsto the time-averaged position of the other of the two beads (404). Thepositions of each bead are projected onto a plane orthogonal to thevector, and the data is kept in 2D-space (406). The two projected beadpositions are spatially averaged with one another, providing for lateralmuscle position in 2D-space (408). The position from 408 is stored, anda linear regression is used to fit a line to the previous N lateralmuscle positions (410). The previous N lateral muscle positions areprojected onto the regression line, and the data is left inone-dimension (412). The one-dimensional data represents the primarydimension of lateral muscle vibration and is stored (414). A FastFourier Transform (FFT) is used to compute the mean peak frequency andamplitude of the previous N primary lateral muscle vibrations (416).

One or more accelerometers (and, optionally, an additional sensor suchas an inertial measurement unit), can be used to identify any vibrationsof the magnetic field sensors relative to the magnetic beads, so as toremove noise from the sensor vibrations from the signal of the lateralvibrations of the one or more magnetic beads.

Magnetic beads can be placed at muscles or tendons to sense lateralmuscle vibrations.

1.4 Shear Wave Tensiometry Via Magnetic Beads Implanted in Muscles,Tendons, or Surrounding Tissue

Methods are provided for measuring a force exerted by or uponmuscle-tendon units. Magnetic targets can be affixed external to orimplanted within tendons for force measurement. Tendon strain measuredvia magnetomicrometry can be employed using a tendon elasticity model toestimate tendon force directly.

A velocity of a shear wave introduced to a tendon can be measured viathe velocity of the shear wave propagation through that tendon. Suchmeasurements are typically obtained via ultrasound measurements at twolocations within the tendon or via measurement of the

acceleration of the tissue superficial to the tendon using twoaccelerometers [26]. Magnetic targets can provide for improvedimplementations of shear wave tensiometry in several regards, using oneor more magnetic beads implanted at the tissue in or around themuscle-tendon unit. One benefit of using magnetic beads instead of usingultrasound or a pair of accelerometers is that magnetic beads can bothprovide high-resolution information in three-dimensions while doing sowith a very small time delay. For example, using one or more magneticbeads, a shear wave can be measured in both transverse dimensions of thetendon, even if they are out of phase with one another. In a furtherexample, two or more magnetic beads can be used, with the time delaybetween the vibrations of the beads being used to determine the velocityof shear waves or compression waves.

A method of estimating a physiological parameter of a muscle or tendonincludes, with a magnetic field sensor (e.g., sensors 106 a-e, or sensorarrays 208 a-c), detecting vibration of a target implanted at a muscleor tendon (e.g., targets 101, 102, 201, 202, 203). The method furtherincludes estimating at least one of muscle force and tendon force basedon the detected vibration.

The method can further include applying a perturbance to the target orto a muscle-tendon unit comprising the muscle and the tendon. A timingof the vibration of the target relative to a timing of the perturbanceor relative to a timing of a vibration of one or more additional targetsimplanted at the muscle-tendon unit can be measured. A speed of a shearwave or a compression wave in the muscle-tendon unit or surroundingtissue of the muscle tendon unit can be estimated based on the measuredtiming. The perturbance can include poking or vibrating themuscle-tendon unit or the surrounding tissue. For example, an appliedmagnetic field, such as can be supplied by an electromagnet, can actuatea magnetic bead affixed at the muscle-tendon unit or the surroundingtissue to initiate the perturbance. A physiological property of themuscle-tendon unit based on the estimated speed of the shear wave or thecompression wave can be determined. The physiological property can be,for example, stiffness of the muscle, the tendon, or the surroundingtissue.

For example, pairs of magnetic beads can be positioned both in thetendon, both in the muscle, one in tendon and one in muscle, or both intissue surrounding or near the tendon or muscle, including on thesurface of the skin external to the body. These waves can be created

using a piezoelectric tapper as in [26], or they can be generated bysome other vibration instrument. Alternatively, the vibrations can becaused by movement of the body or by the lateral vibrations of themuscle during activation, and, optionally, this muscle movement can besensed using electromyography or magnetic-bead-sensed lateralvibrations, as described above. In an example method, a time delaybetween the source of the compression or shear wave (whether from atapping device, muscle movement, electromagnet, or some other source)and the measurement of the vibration (whether via a magnet in themuscle, tendon, or surrounding tissue, including on the surface of theskin external to the body) can be used to determine the velocity of thewave. In a further example method, magnetic beads can be used todistinguish between mechanomyographic lateral vibrations, shear waves,and compression waves. For example, because mechanomyographic lateralvibrations are associated with muscle fibers and shear waves areassociated with tendon fibers, the differences in physiology can giverise to differences in parameters such as wave velocity or wavelength,allowing the two measurements to be performed using overlapping sets ofsensors. An electromagnetic coil can be used to vibrate a magnetic beadto create the shear wave or compression wave.

Another advantage provided by the use of magnetic beads is that thereflection from a transient wave, albeit smaller than an original wave,may be sensed with greater certainty. This allows for an ability tomeasure wave velocity using a single magnetic bead, or to provide arefined measurement when measuring wave velocity using multiple beads.Or in a further example, this can allow the use of a magnetic bead nearthe musculotendinous junction to be used to measure delay between thetraveling of a tranverse wave from the muscle into the tendon and itsreflection back into the muscle after traveling and returning the lengthof the tendon.

1.5 Applications

The use of magnetic targets, as described above, can provide forimprovements and new capabilities in biomechatronics, including, forexample improving prosthetic control, both in controlling robotic hands,as well as in controlling robotic legs and feet. Such methods anddevices can also provide for improved control over orthotic andexoskeletal devices for support, rehabilitation, and augmentation.Additional applications include untethered animal and human biomechanicsstudies in natural environments, which can elucidate new principles ofmotion; virtual reality hand and tool tracking, which can benefit fromthis new tracking strategy due to its ability to provide fast, accurateposition information without the need for line-of-sight; and, sensing ofmuscle lengths and velocities in real time, which can enable closed-loopmuscle stimulation control feedback, opening the door to improvedstrategies for paralysis mitigation and recovery.

1.5.1 Control of Robotics Via Magnetic Beads

There are various ways that implanted magnetic beads can be used for thecontrol of robotics, such as wearable robotics such as prosthesis,orthoses, and exoskeletons, or other external robotic devices such ashumanoids, cars or power tools. An example control diagram is providedshowing how magnetomicrometry and mechanomyography can be used in thecontrol of, specifically, a robotic prosthesis (see FIG. 17 ). However,the same or similar control methodologies can be used for the control ofan exoskeleton, orthosis or other such device.

As shown in FIG. 17 , tracking of magnetic beads (502), as describedabove, can provide for muscle length measurements via magnetomicrometry(504) and muscle activation measurements via mechanomyography (508).With the monitoring of time (506), muscle velocity can also be obtained(512). Any or all of the obtained muscle length, muscle velocity, andmuscle activation can be provided to a biophysical muscle model (e.g.,the Hill Muscle Model) capable of predicting muscle force given inputsof muscle length, velocity and activation (510). The obtained muscleforce can then be provided for force control of a robotic prosthetists(516). Force and/or position measurements obtained from the prosthesiscan provide for closed-loop control (518) by providing feedback to abiophysical model of the anatomy (514) and to the force control paradigm(516).

Of course, there are many variations with which this control diagram canbe used. For example, instead of using force control (516) to control arobotic prosthesis (e.g., prosthesis 210), the muscle length (504) orvelocity (512) can be used to control the device using position orvelocity control, respectively, in a master-slave configuration, or theforce can be used to control the position or the velocity of the device.Alternatively, or in addition, electromyography can provide formechanomyography in this control scheme. Any of these control methodscan be performed with either open-loop or closed-loop control (e.g.,feeding back information from the prosthesis as to position, velocity,or force), or using control methodologies, such as Kalman filtering, inthe control. Any subset or combination of the inputs, outputs, andmethodologies shown in FIG. 17 can be provided for control of a roboticprosthesis.

Furthermore, biological feedback mechanisms such as electricalstimulation, magnetic bead vibration via external coils, or tactilefeedback via skin pressure or vibration, can be included in the controlof the device. Further still, a measured length, velocity, and force ofthe muscle can be combined with a biophysical model of the limbincluding inertia terms to determine biomechanical information, such asjoint states and torques, about the user of the device, which can thenbe fed back into the control algorithm for closed-loop kinetic orkinematic control.

1.5.2 Portable Biomechanics Monitoring Using Magnetic Beads—PortableGoniometry

Magnetic targets can be used for motion capture biomechanical studies“in the wild.” Due to the traditional need for an array of cameras, eachof which requires line-of-sight to biomechanical markers, motion capturedata collections are typically confined to a small lab space. Thefreedom to have a biomechanical subject (human or animal) move about ina natural setting imparts the ability to study the subject in a morerelevant environment and also allows for the biomechanical signals to beused for additional human-computer interfacing. Further, the ability tomonitor biomechanical motion capture markers without a need for line ofsight allows for the subject to wear comfortable, loose fitting clothingduring biomechanical monitoring, encouraging the subject to move in amore natural manner.

An example of a magnetic bead position sensing system and method thatcan provide for freedom in biomechanical motion capture monitoring isshown in FIG. 18 . As shown in FIG. 18, the system 600 includes magnetictargets 601, 602 disposed at a joint (e.g., the hip) of a subject. Amagnetic field sensing array 608 is disposed nearby the targets.

A method of monitoring biomechanical motion includes disposing at leastone target (601, 602) on a subject at a location associated with a jointof the subject, and, with a magnetic field sensor array (608), detectinga change in state of the at least one target relative to the magneticfield sensor array, relative to another target disposed on the subject,or a combination thereof. A state of the joint can then be determinedbased on the detected change in state of the at least one target.

The human body has a diverse array of joints, with differing degrees offreedom and differing ranges of motion, thus requiring an array ofdifferent sensing options. Systems and methods involving magnetic targettracking can take many different forms to account for this diversity inability of movement. It is not sufficient, for instance, to use aprotractor to describe the position or the orientation of the shoulderand upper arm, because the shoulder allows the arm to move up and down(abduction/adduction), as well as forward and back (flexion/extension),as well as in rotation.

At least one magnetic bead can be used as a biomechanical marker totrack a position or orientation of a joint for biomechanical monitoring.When the magnetic field sensing array is used as a biomechanical markeritself, the one or more beads can be tracked in position and orientationrelative to the sensing array and to each other. When the magnetic fieldsensing array is not itself used as a biomechanical marker, the two ormore beads can be tracked relative to one another. The biomechanicalmarkers (e.g., the sensing array, if used as a marker, as well as theone or more magnetic beads) can be affixed at a body segment. Forexample, the markers can be fixed using an adhesive, such as tape, orthe markers can be manufactured as a component of a sticker whichaffixes to the surface of the skin, or the markers can be embedded inclothing (e.g., tightly worn clothing). To enable the biomechanicalmonitoring to be performed more accurately or with fewer sensors, theorientations of the magnetic markers can be indicated to ensure optimalplacement. For every position or orientation of the markers relative toone another, the position and orientation of the body segments relativeto one another can be calculated.

As one example (see FIG. 18 ), the human hip can be monitored with twodegrees of freedom using a magnetic marker placed directly on thesurface of the skin over the lateral edge of the ilium (the outside ofthe hip bone), with the north pole directed anteriorly (directedforward, or toward the front of the body), a magnetic marker placeddirectly on the surface of the skin laterally to the femur toward theproximal end of the upper leg, with the north pole directed distally(directed down the leg, or toward the knee), and a magnetic fieldsensing array not used as a biomechanical marker (so that it is free to“float”) can be placed outside of clothing between the two markers tomonitor their positions relative to one another. The relative positionsand angles of the two magnetic markers can then be used to identify theposition of the femur relative to the pelvis.

Methods and systems for magnetomicrometry, as described herein, canprovide for tracking of the target(s) and/or sensor array(s). Fordetermination of a joint state, tracking information obtained can becombined with information relating to the anatomical structure of theportion of the body at which the targets are disposed, including, forexample, the type of joint and a number of degrees of freedom of thejoint.

1.6 Magnetic Bead Insertion Device

As discussed in Example 1, in this early feasibility study, a needle anda pair of surgical scissors were used to make an insertion channel forthe magnet, and a hollow plastic tube was used to insert the magnet withthe aid of a wooden push rod, after which the muscle and skin weresutured closed and film dressing was applied. In clinical practice,there are instances in which the magnets may be inserted during asurgical procedure with the muscle similarly exposed. To generalize thisprocedure, a device is provided which is capable of additionallyinserting the magnetic beads through the skin boundary without the needof a surgical procedure that exposes the muscle. The device can allowfor either percutaneous or surgical implantation of one or more magneticbeads.

An insertion device 700, alternatively referred to as an injector, isshown in FIG. 19 . The device includes a cannula 705 and a cartridge 707configured to position an implantable target 701 at the cannula. Apushrod 711 is received in the cannula and configured to push theimplantable target 701 from the cartridge through the cannula fordelivery to a tissue. The cannula and the pushrod comprise a nonmagneticmaterial. A distal end 709 of the cannula can be of a complementarygeometry to the implantable target to prevent damage to a shellstructure of the target during delivery. For example, where the targetis a spherical magnetic bead, the distal end of the cannula can be of acurved geometry. A distal surface 712 of the pushrod can include aspring structure 714 that can prevent damage to the shell structure ofthe target (e.g., an elastomeric material, a coil spring, etc.). Thedevice can further include a mount 715 coupling the cartridge 707 to thecannula. The cartridge can be configured to retain a plurality oftargets and can be rotatable about the mount to position each of theplurality of targets at the cannula.

An insertion kit or an insertion system can include a plurality ofimplantable targets (e.g., magnetic beads 300) and an insertion device(e.g., device 700). The implantable targets can be deliverable by theinsertion device.

As described above in the discussion of fluoromicrometry, nonmagneticbeads have long been used for tracking tissue positions using X-rays, atechnique called radiosteriometric analysis (RSA). To perform RSA, thenonmagnetic beads (typically made of tantalum) are injected, oftenpercutaneously, into the tissue to be tracked. This injection process isperformed using devices such as the tantalum bead inserter manufacturedand sold by RSA Biomedical and Halifax Biomedical.

Magnetic bead tracking differs in that the injected beads aremagnetized, and that they are often of larger diameter. Devices, such asinjector 700, can inject permanently magnetized implantable components,such as spherical, cylindrical, or cube magnets (or any other geometryof magnet), into human or animal tissue via a nonmagnetic needle (e.g.,barrel, cannula, or rigid tube) and using a nonmagnetic pushrod. Theneedle and the pushrod can be of a same or different material, composedof rigid nonmagnetic materials, such as titanium, copper, nonmagneticstainless steel, gold, wood, silver, brass, glass, aluminum, zinc,marble, bronze, or a polymer (such as a plastic). The term “nonmagnetic”is defined herein to mean any material that is not temporarily orpermanently magnetizable, such that the material does notelectromagnetically interact with a permanently magnetized object whenthe material and the object are static relative to one another and noelectric current is introduced into the system. The pushrod or needlecan be coated in a biocompatible material to increase thebiocompatibility of the surface and to decrease the stiffness of thesurface, to decrease the possibility of the biocompatible coating of theinserted component being compromised (for instance, with a coating ofgold or parylene).

One or more implantable magnetic components can be inserted into anonmagnetic cartridge, such as in the revolver or magazine of a firearm,wherein the cartridge is inserted into the device for the injection ofthe one or more magnetic implants. The cartridge can have magnetized orferromagnetic components to guide the magnetic flux, to help preservethe magnetization of the permanently magnetized components over time(days to decades), or to hold the beads in a particular location ororientation before being dislodged by the pushrod. The pushrod caninclude a spring in order to minimize the possibility of excessive forcebeing placed on the bead. The pushrod can have a smooth tip that matchesthe geometry of the implant being injected into the body, or that isflat, beveled, or domelike. The above instances can be adapted toaccommodate magnetic implants of various sizes, for instance sphericalmagnetic beads of 1 mm, 2 mm, 3 mm, 4 mm, or 5 mm diameter, orcylindrical or cube magnets in a similar range of sizes.

As an example of several of these instances being combined, FIG. 19shows a plastic pushrod with a curved tip used to dislodge a sphericalmagnetic bead from a cartridge of multiple beads and guide it down anonmagnetic stainless-steel cannula into a muscle, though an insertionpath that has been created using the sharp end of the cannula. In FIG.19 , the cartridge is mounted so that it can revolve with respect to thecannula, and the pushrod is not rigidly affixed to the other componentsof the device.

The magnetic bead insertion device can be used in combination with animaging device, such as ultrasound, an augmented reality device combinedwith previously-collected Mill data, or a magnetic bead tracking system,to guide the insertion of the magnetic beads.

2. Multi-Valued Inverse Existence Proof for N Sensors 2.1 Proof ofNon-Guarantee of Single-Valued Inverse

It is shown that the inverse of the magnetic dipole field equations isnot guaranteed to be single-valued. To do this, a demonstration isprovided of the existence of at least one geometry, for any number ofsensors and any magnet strength, where the magnet can cause identicalmagnetic fields at all of the sensors from multiple, distinct locations.

Reference is made herein to the Taylor Reference [17], the entireteachings of which are incorporated herein by reference.

Consider N sensors on a circle, not necessarily evenly spaced, and amagnet oriented normal to the circle, positioned on a line normal to andintersecting the center of the circle, as illustrated in FIG. 20 .

To investigate this situation, we first convert equations (8a) and (8c)of the Taylor Reference into cylindrical coordinates, and with y⁻=0 andθ=φ=0, we have

$\begin{matrix}{B_{\rho} = {\overset{\_}{m}\left( \frac{3\overset{\_}{\rho}\overset{\_}{z}}{\left( {{\overset{\_}{\rho}}^{2} + {\overset{\_}{z}}^{2}} \right)^{5/2}} \right)}} & \left( {2.1a} \right)\end{matrix}$ $\begin{matrix}{{B_{z} = {\overset{\_}{m}\left( {\frac{3{\overset{\_}{z}}^{2}}{\left( {{\overset{\_}{\rho}}^{2} + {\overset{\_}{z}}^{2}} \right)^{5/2}} - \frac{1}{\left( {{\overset{\_}{\rho}}^{2} + {\overset{\_}{z}}^{2}} \right)^{3/2}}} \right)}},} & \left( {2.1b} \right)\end{matrix}$

where ρ⁻ is a vector giving any sensor's radial position relative to themagnet, and B_(ρ) is the radially-directed magnetic field sensed at allof the sensors. For simplicity, we consider the case where the magneticfield from the magnet is fully in the plane of the sensors at each ofthe sensors (i.e., B_(z)=0). We then have

$\begin{matrix}{0 = {\overset{\_}{m}\left( {\frac{3{\overset{\_}{z}}^{2}}{\left( {{\overset{\_}{\rho}}^{2} + {\overset{\_}{z}}^{2}} \right)^{5/2}} - \frac{1}{\left( {{\overset{\_}{\rho}}^{2} + {\overset{\_}{z}}^{2}} \right)^{3/2}}} \right)}} & (2.2)\end{matrix}$ $\begin{matrix}{= {\overset{\_}{m}\left( \frac{{3{\overset{\_}{z}}^{2}} - \left( {{\overset{\_}{\rho}}^{2} + {\overset{\_}{z}}^{2}} \right)}{\left( {{\overset{\_}{\rho}}^{2} + {\overset{\_}{z}}^{2}} \right)^{5/2}} \right)}} & (2.3)\end{matrix}$ $\begin{matrix}{= {\overset{\_}{m}\left( \frac{{2{\overset{\_}{z}}^{2}} - {\overset{\_}{\rho}}^{2}}{\left( {{\overset{\_}{\rho}}^{2} + {\overset{\_}{z}}^{2}} \right)^{5/2}} \right)}} & (2.4)\end{matrix}$ $\begin{matrix}{\left. \Rightarrow 0 \right. = {{2{\overset{\_}{z}}^{2}} - {\overset{\_}{\rho}}^{2}}} & (2.5)\end{matrix}$ $\begin{matrix}{\overset{\_}{z} = {{\pm \frac{1}{\sqrt{2}}}{\overset{\_}{\rho}.}}} & (2.6)\end{matrix}$

Thus, when the magnet is positioned above or below the circle of sensorsat a height or depth of approximately 0.707 times the radius of thecircle of sensors, the z-directed magnetic field seen by the sensors iszero. Because the magnet is oriented with the north pole facing up, ifthe magnet is positioned above the array, the magnetic fields aredirected radially inwards, but if the magnet is positioned below thearray, the magnetic fields are directed radially outwards with the samemagnitude. However, reversing the polarity of the magnet reverses theradial direction of the magnetic field without changing the magnitude.For example, a magnet facing north-pole-up at a height above the sensorsof 0.707 times the sensor circle radius will produce the same magneticfields at the sensors as a magnet facing north-pole-down at a depth of0.707 times the radius under the sensors (see FIG. 20 ).

Described mathematically, this is shown simply as

$\begin{matrix}{B_{\rho} = {\overset{\_}{m}\left( \frac{3\overset{\_}{\rho}\overset{\_}{z}}{\left( {{\overset{\_}{\rho}}^{2} + {\overset{\_}{z}}^{2}} \right)^{5/2}} \right)}} & (2.7)\end{matrix}$ $\begin{matrix}{= {\overset{\_}{m}\left( \frac{3{\overset{\_}{\rho}\left( {\frac{1}{\sqrt{2}}\overset{\_}{\rho}} \right)}}{\left( {{\overset{\_}{\rho}}^{2} + \left( {\frac{1}{\sqrt{2}}\overset{\_}{\rho}} \right)^{2}} \right)^{5/2}} \right)}} & (2.8)\end{matrix}$ $\begin{matrix}{= {\left( {- \overset{\_}{m}} \right)\left( \frac{3{\overset{\_}{\rho}\left( {{- \frac{1}{\sqrt{2}}}\overset{\_}{\rho}} \right)}}{\left( {{\overset{\_}{\rho}}^{2} + \left( {{- \frac{1}{\sqrt{2}}}\overset{\_}{\rho}} \right)^{2}} \right)^{5/2}} \right)}} & (2.9)\end{matrix}$ $\begin{matrix}{= {{\overset{\_}{m}}^{\prime}\left( \frac{3\overset{\_}{\rho}{\overset{\_}{z}}^{\prime}}{\left( {{\overset{\_}{\rho}}^{2} + {\overset{\_}{z}}^{\prime 2}} \right)^{5/2}} \right)}} & (2.1)\end{matrix}$

where m bar prime and z bar prime represent, respectively, the magneticdipole weight of the magnet and the vertical position of the sensorsrelative to the magnet.

Having shown that, for this particular geometry, at least onemulti-valued inverse exists for any strength magnet regardless of thenumber of sensors used to track the magnet, it is thus concluded thatthe inverse is not guaranteed to be single-valued for all geometries,regardless of the number of sensors used.

Additional configurations in which multiple, distinct magnet positionsyield the same magnetic field at multiple sensor positions are possible.For instance, when tracking magnets with six degrees of freedom in thecircular sensor configuration discussed above, for each different magnetheight, a phantom magnet of differing magnetic dipole weight with adepth that differs from the magnet height can appear. As a furtherexample, any locations relative to the magnet where the magnetic fieldvectors are both colinear and either parallel or anti-parallel cancreate a rotational degree of freedom about a pair of sensors when thecolinear parallel or anti-parallel magnetic fields correspond with thesensor locations.

Now, perhaps it can be proved that a single-valued inverse is guaranteedwhen more than two sensors are used that are not at coplanar locationson a circle. However, the existence proof above is still relevant forthese particular cases because it reveals locations of potential localminima and suggests how sensor geometries can be designed to avoid theseminima. It also suggests the difficulty in finding an analytic inverseto the magnetic dipole field equation, as the analytic inverse wouldneed to take into account sensor geometry and the possible presence ofmultiple solutions.

The use of the spatial magnetic field gradient has been shown to allowinversion from magnetic dipole field measurements, giving dipolelocations via matrix inversion [23, 10, 22], but these direct-inversionmethods require calculating the spatial gradient of the field, which,with presently-available magnetic field sensors, introducesnon-negligible numerical error. More importantly, the direct-inversiontechniques are currently limited to just a single magnet tracked via asingle estimated magnetic field and gradient, calling for furtherrefinement via optimization methods. Extending these direct-inversiontechniques to the tracking of multiple magnets using many magnetic fieldsensors can improve the accuracy of direct-inversion using the spatialgradient and perhaps providing a high-speed, high-accuracy global searchstrategy for many magnets at once. Modifications to compensate formagnetic disturbance fields can be provided.

3. Measurement of Dipoles

The following section presents improvements to the methods and systemsrelating to tracking of magnetic objects found in WO2019/074950, “Methodfor Neuromechanical and Neuroelectromagnetic Mitigation of LimbPathology” and in the Taylor Reference, the entire teachings of whichare incorporated herein by reference.

A method for determining one or more of three sensor position parametersand three sensor orientation parameters for each of the sensors in asensor array includes placing at least one target in at least one knownlocation relative to a sensor array, whereby a signal from the at leastone target at the sensors is detected, and recording at least onemeasurement of the signal at each of the sensors for each placement ofthe one or more targets. The method further includes estimating one ormore parameters from the group consisting of x-position, y-position,z-position, yaw, pitch, and roll, of each of the sensors and estimatingany unknown state parameters of the at least one target. A constantvalue for a magnetic dipole weight state parameter of the at least onetarget for each of the measurements is provided, and predicted values ofthe signal at each of the sensors for each of the measurements given theone or more estimated sensor parameters, the estimated target stateparameters, and the provided constant magnetic dipole weight stateparameter are calculated. A prediction error in the predicted values ofthe signal with reference to the values of the signals detected at thesensors is computed. A prediction error Jacobian matrix by analyticallycomputing elements of the prediction error Jacobian matrix with respectto the estimated parameters of the sensors for each measurement iscalculated and, from the prediction error and the prediction errorJacobian matrix, a state of the parameters of the sensors is determined.

3.1 Strength Initialization

In the Taylor Reference, Eqn. 7 was used to calculate an initialestimate of the magnetic dipole weight, followed by a short (˜30 s)six-degree-of-freedom tracking session in which the magnetic dipoleweight was continuously tracked and recorded. The median of the recordedmagnetic dipole weights for each magnet was then used as the knownmagnetic dipole weight value for that magnet in five-degree-of-freedomtracking. For maximum stability, accuracy, and bandwidth, each magnet'sstrength should be measured before tracking (as was performed in thiswork) whenever an application permits doing so.

It was assumed in the Taylor Reference that the residual flux density(the remanence), B_(rj), of a magnet is readily accessible information.This parameter is not, however, always reported by magnet vendors.Strategies are shared for using either the “N-Rating” or the surfacefield of a magnet to determine an initial estimate of the magnet'sresidual flux density.

In practice, in this work it was found to be tolerable for an initialestimate of the dipole strength to be incorrect by several orders ofmagnitude when starting six-degree-of-freedom tracking (i.e., whenmeasuring the magnet's strength). However, in multiple-magnetsix-degree-of-freedom tracking (e.g., when using a sensor array tomeasure the strength of two magnets at once), a more accurate initialestimate can improve stability at the start of tracking. Additionally,the formulas below aid in performing simulations to determine whichfactory-specified magnets will work best for a given application.

For optimal accuracy when sensor repeatability is in doubt, the magneticdipole weight can be measured using the same sensors and calibrationthat can be used in five-degree-of-freedom tracking, since the strengththat is measured by a sensor array is measured relative to thesensitivity of the given sensors.

Although the methods below are derived for spherical magnets, these samestrategies can be used to determine an approximation for the strength ofa magnet of similar geometry, such as a cube or a cylinder, forfar-field tracking.

3.1.1 Strength from N-Rating

If the residual flux density is known, an initial guess for the magneticdipole weight can be calculated via Eqn. 7 in the Taylor Reference. Ifeither the N-rating or the magnetization of the magnet are known, thesecan be used to estimate the residual flux density. (If the magnetizationM_(j) is known

$\left( {{{where}M_{j}} = \frac{m_{j}}{V_{j}}} \right),$

the residual flux density can be calculated from it as B_(rj)=μ₀M_(j).)Using tables from several prominent magnet vendors [6, 7, 19], the ruleof thumb can be empirically-derived, where N_(j) is the N-rating of thejth magnet.

$\begin{matrix}{B_{rj} = \frac{\frac{N_{j}}{6} + 6}{10}} & (3.1)\end{matrix}$

However, if the manufacturer is known and the manufacturer has a tablereporting the residual flux density, the manufacturer's residual fluxdensity specifications can be used instead. The magnet specifications asreported by a manufacturer may not line up with the dipole magnitude asseen by the sensors, due to errors in the magnet specifications, thesensor specifications, or in the above empirically-derived formula.

3.1.2 Strength from Surface Field

Alternatively, an estimate of the strength of the magnet can bedetermined from the surface field of the magnet measured at its north orsouth pole by a gaussmeter (also known as a fluxmeter or teslameter, agaussmeter is a magnetometer with a very high full-scale range) or asreported by the factory specifications of a given magnet. The residualflux density of a spherical magnet is not the same as its field at thenorth pole. Repeating Eqn. 8c of the Taylor Reference here forconvenience, the z-directed field is given by:

$B_{iz} = {G_{z} + {\sum\limits_{j^{\prime} = 1}^{j^{\prime} = M}{{\overset{\_}{m}}_{j^{\prime}}\left( {\frac{3{{\overset{\_}{z}}_{{ij}^{\prime}}\left( {{\sin\theta_{j^{\prime}}\cos\phi_{j^{\prime}}{\overset{\_}{x}}_{{ij}^{\prime}}} + {\sin\theta_{j^{\prime}}\sin\phi_{j^{\prime}}{\overset{\_}{y}}_{{ij}^{\prime}}} + {\cos\theta_{j^{\prime}}{\overset{\_}{z}}_{{ij}^{\prime}}}} \right)}}{\left( {{\overset{\_}{x}}_{{ij}^{\prime}}^{2} + {\overset{\_}{y}}_{{ij}^{\prime}}^{2} + {\overset{\_}{z}}_{{ij}^{\prime}}^{2}} \right)^{5/2}} -} \right.}}}$$\left. \frac{\cos\theta_{j^{\prime}}}{\left( {{\overset{\_}{x}}_{{ij}^{\prime}}^{2} + {\overset{\_}{y}}_{{ij}^{\prime}}^{2} + {\overset{\_}{z}}_{{ij}^{\prime}}^{2}} \right)^{3/2}} \right)$

The surface field at the north pole of the jth spherical magnet iscalculated by setting θ_(j)=⁻x_(ij)=⁻y_(ij)=G_(z)=0 and z⁻ _(ij)=R_(j)in this equation, giving

$\begin{matrix}\begin{matrix}{B_{sj} = {{\overset{\_}{m}}_{j}\left( {\frac{3{{\overset{\_}{z}}_{ij}\left( {{\sin\theta_{j}\cos\phi_{j}{\overset{\_}{x}}_{ij}} + {\sin\theta_{j}\sin\phi_{j}{\overset{\_}{y}}_{ij}} + {\cos\theta_{j}{\overset{\_}{z}}_{ij}}} \right)}}{\left( {{\overset{\_}{x}}_{ij}^{2} + {\overset{\_}{y}}_{ij}^{2} + {\overset{\_}{z}}_{ij}^{2}} \right)^{5/2}} -} \right.}} \\\left. {}\frac{\cos\theta_{j}}{\left( {{\overset{\_}{x}}_{ij}^{2} + {\overset{\_}{y}}_{ij}^{2} + {\overset{\_}{z}}_{ij}^{2}} \right)^{3/2}} \right) \\{= {{\overset{\_}{m}}_{j}\left( {\frac{3{R_{j}\left( R_{j} \right)}}{R_{j}^{5}} - \frac{1}{R_{j}^{3}}} \right)}} \\{= {{\overset{\_}{m}}_{j}\frac{2}{R_{j}^{3}}}}\end{matrix} & (3.2)\end{matrix}$

Then, comparing (3.2) with Eqn. 7 of the Taylor Reference gives:

B_(rj)= 3/2B_(rj)  (3.3)

3.1.3 Other Methods

There are other methods to determine the strength of a magnet (forinstance, the time integral of a coil rotated about the magnet can beused to determine the magnet strength). Moreover, if the magnet isnon-spherical, the above methods can still be used to determine anacceptable initialization. The estimation does not have to be precise inorder to begin six-degree-of-freedom tracking. For a single magnet, theinitial estimate can be incorrect by several orders of magnitude andstill be tracked. However, for multiple magnets, it becomes morenecessary, for stability, to have better initial estimates.

3.2 Non-Negative Strength Constraint in Six-Degree-of-Freedom Tracking

On a practical note, it was found in this work that the trackingalgorithm was more stable if negative magnetic dipole strengths werereversed in sign and orientation during six-degree-of-freedom tracking,especially when tracking multiple magnets at once. In other words, itwas found that tracking stability was increased by multiplying m⁻ _(j)by −1 and subtracting π from θ_(j) whenever the tracking algorithmdetermined that m⁻ _(j)<0. This was found to be particularly importantin the stability of multi-magnet measurement.

3.3 Refinement to Dipole Measurement

It is possible to refine the magnetic dipole weight measurement beyondmerely using the median of the recorded magnetic dipole weights duringsix-degree-of-freedom tracking. The measurement process can first beperformed as described above and in view of the Taylor Reference, but,in addition to recording the magnetic dipole weight, the tracking data(e.g., position and orientation of the magnets) can be recorded as well.The median of the magnetic dipole weight, along with thefive-degree-of-freedom data across all (or a representative subset) ofthe tracking timesteps can then be used as an initial estimate for a newoptimization. This new optimization can be performed similarly tocurrent tracking methods, except that instead of the error correspondingto a single timestep over all sensors, the optimization takes intoaccount the error across all timesteps over all sensors. In other words,both the path and the magnetic dipole weight can be optimized, and theoutput can be a refined path with a single magnetic dipole weightcorresponding to each magnet.

A description of example optimization is provided. As before, themagnetic field prediction error is used for the optimization costfunction. At the ith sensor and kth timestep, the magnetic fieldprediction error, E_(jk), is the difference between the measuredmagnetic field B^({tilde over ( )}) _(jk) and the predicted magneticfield B_(jk),

E _(jk) =B _(jk) −{tilde over (B)} _(jk),  (3.4)

For each timestep (the kth timestep), the derivatives of the magneticfield prediction errors with respect to each of the full-path magneticdipole weight estimates are calculated and placed in a separatesubmatrix

$\begin{matrix}{M_{k} = {\begin{bmatrix}{\frac{\partial}{\partial{\overset{\_}{m}}_{1}}B_{1{kx}}} & & {\frac{\partial}{\partial{\overset{\_}{m}}_{M}}B_{1{kx}}} \\{\frac{\partial}{\partial{\overset{\_}{m}}_{1}}B_{1{ky}}} & \ldots & {\frac{\partial}{\partial{\overset{\_}{m}}_{M}}B_{1{ky}}} \\{\frac{\partial}{\partial{\overset{\_}{m}}_{1}}B_{1{kz}}} & & {\frac{\partial}{\partial{\overset{\_}{m}}_{M}}B_{1{kz}}} \\{\frac{\partial}{\partial{\overset{\_}{m}}_{1}}B_{2{kx}}} & & {\frac{\partial}{\partial{\overset{\_}{m}}_{M}}B_{2{kx}}} \\{\frac{\partial}{\partial{\overset{\_}{m}}_{1}}B_{2{ky}}} & \ldots & {\frac{\partial}{\partial{\overset{\_}{m}}_{M}}B_{2{ky}}} \\{\frac{\partial}{\partial{\overset{\_}{m}}_{1}}B_{2{kz}}} & & {\frac{\partial}{\partial{\overset{\_}{m}}_{M}}B_{2{kz}}} \\ \vdots & \ddots & \vdots \\{\frac{\partial}{\partial{\overset{\_}{m}}_{1}}B_{Nkx}} & & {\frac{\partial}{\partial{\overset{\_}{m}}_{M}}B_{Nkx}} \\{\frac{\partial}{\partial{\overset{\_}{m}}_{1}}B_{Nky}} & \ldots & {\frac{\partial}{\partial{\overset{\_}{m}}_{M}}B_{Nky}} \\{\frac{\partial}{\partial{\overset{\_}{m}}_{1}}B_{Nkz}} & & {\frac{\partial}{\partial{\overset{\_}{m}}_{M}}B_{Nkz}}\end{bmatrix}.}} & (3.5)\end{matrix}$

Defining J_(k) to be the Jacobian matrix corresponding to the kthmeasurement, where the elements corresponding to the full-path magneticdipole weight have been removed and placed in M_(k), the full-pathdipole measurement Jacobian matrix is then constructed as

$\begin{matrix}{\overset{\_}{J} = {\begin{bmatrix}J_{1} & 0 & \ldots & 0 & M_{1} \\0 & J_{2} & \ldots & 0 & M_{2} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\0 & 0 & \ldots & J_{K} & M_{K}\end{bmatrix}.}} & (3.6)\end{matrix}$

For example, with a single magnet measured by N sensors over Ktimesteps, this is constructed as shown in FIG. 35 , where x⁻ _(ijk), y⁻_(ijk), and z⁻ _(ijk) are the distances in x, y, and z to the ith sensorfrom the jth magnet at the kth timestep, and θ_(jk) and φ_(jk) are theorientation from and around vertical of the jth magnet at the kthtimestep.

Note that, as with the magnetic dipole weight in the example above, anyother tracking parameter can be considered constant (with either anunknown or a previously known value) across some or all measurements.For instance, the position can be fixed while the orientation is allowedto change (e.g., the magnet being rotated while positioned in anindentation), the orientation can be fixed while the position is allowedto change (e.g., a magnet fixed in a block that is moved across asurface), and disturbance field of any given complexity can be accountedfor and considered varying or constant across the measurements.

FIG. 21 is a flowchart illustrating an example method of performing arefined dipole measurement. The method 800 includes beginning with anestimate of the strength of the magnet and its initial location andorientation (802). For K timesteps, at each timestep (804): the magneticfield at the N sensors is measured (806), the parameters of the magnet(e.g., location, orientation, and strength) are estimated (808), themagnetic field measurements and magnet parameter estimates are recorded(810), and the magnet is translated or rotated (812). Using all orseveral recorded magnetic field measurements and magnet parameterestimates, the locations and orientations of the magnet across all Ktimesteps are re-estimated while estimating a single magnet strengthparameter (e.g., a constant value) (814). The single magnet strengthparameter is recorded for use in subsequent magnet tracking (816).

3.4 Minimum Distance of a Magnet from any Magnetometer

It can be prudent in considering the measurement and tracking of amagnet to also consider the full-scale range of the sensors that areused to measure and track the magnet. In this section, strategies areprovided for making a conservative estimate of the minimum distance thatis to be maintained between a single permanent magnet and any magneticfield sensor while not exceeding the full-scale range of the sensor.

To calculate a conservative estimate, the equation for magnetic field isreduced to a function of distance from the magnet along the axis thathas the greatest magnetic field. The scenario is then considered inwhich this axis of the magnet is aligned with one of the three axes ofthe sensor. The axis of a magnet along which the magnetic field isgreatest at a given radial distance is the z-axis of the magnet, so Eqn.8c of the Taylor Reference is used, with M=N=1 and x⁻=⁻y=θ=G_(z)=0 at adistance of z⁻=d_(min):

$\begin{matrix}{\begin{matrix}{G_{\max} = {G_{z} + {\sum\limits_{j^{\prime} = 1}^{j^{\prime} = M}{\overset{\_}{m}}_{j^{\prime}}}}} \\\left( {\frac{3{{\overset{\_}{z}}_{{ij}^{\prime}}\left( {{\sin\theta_{j^{\prime}}\cos\phi_{j^{\prime}}{\overset{\_}{x}}_{{ij}^{\prime}}} + {\sin\theta_{j^{\prime}}\sin\phi_{j^{\prime}}{\overset{\_}{y}}_{{ij}^{\prime}}} + {\cos\theta_{j^{\prime}}{\overset{\_}{z}}_{{ij}^{\prime}}}} \right)}}{\left( {{\overset{\_}{x}}_{{ij}^{\prime}}^{2} + {\overset{\_}{y}}_{{ij}^{\prime}}^{2} + {\overset{\_}{z}}_{{ij}^{\prime}}^{2}} \right)^{5/2}} -} \right. \\\left. {}\frac{\cos\theta_{j^{\prime}}}{\left( {{\overset{\_}{x}}_{{ij}^{\prime}}^{2} + {\overset{\_}{y}}_{{ij}^{\prime}}^{2} + {\overset{\_}{z}}_{{ij}^{\prime}}^{2}} \right)^{3/2}} \right) \\{= {\frac{B_{r}}{3}{R^{3}\left( {\frac{3d_{\min}^{2}}{d_{\min}^{5}} - \frac{1}{d_{\min}^{3}}} \right)}}} \\{= {\frac{2}{3}\frac{B_{r}R^{3}}{d_{\min}^{3}}}}\end{matrix},} & (3.7)\end{matrix}$

where Eqn. 7 of the Taylor Reference was used to reparameterize themagnetic field in terms of the magnet radius, R. Rearranging (3.7) thengives:

$\begin{matrix}{{d_{\min} = {\sqrt[3]{\frac{2}{3}\frac{B_{r}}{B_{\max}}} \cdot R}},} & (3.8)\end{matrix}$

where d_(min) and R have the same units as one another, and B_(r) andB_(max) have the same units as one another. For a given magnet residualflux density and given sensor, then, the minimum distance that a singlemagnet can be tracked from that sensor is proportional to the magnetradius.

The LSM9DS1 and LIS3MDL magnetic field sensors (STMicroelectronics, NV)have a full-scale range of B_(max)=(2⁽¹⁶⁻¹⁾−1)·58·10⁻⁹ T, or 1.9 mT, sowith these sensors it may not possible to track a magnet at a distancecloser than the d_(min) corresponding to this B_(max). However, theminimum distances corresponding to the recalibration point, B_(max)=5mT, as well as to the maximum exposed field, B_(max)=100 mT (beyondwhich field level permanent damage to the sensors may occur), can alsobe considered. For magnets of various residual flux densities, Table 3.1lists the minimum distances, in multiples of the magnet radius,corresponding to each of the above magnetic field limits when tracking asingle magnet.

TABLE 3.1 Minimum Magnet Distance from a Single Magnet to an LSM9DS1 orLIS3MDL Sensor in Multiples of Magnet Radius, Given Residual FluxDensity Residual Flux Density (T) 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5(N-Rating) (N12) (N18) (N24) (N30) (N36) (N42) (N48) (N54) TrackingMinimum 6.55 6.81 7.06 7.29 7.50 7.70 7.90 8.08 Recalibration Point 4.754.94 5.11 5.28 5.43 5.58 5.72 5.85 Maximum Exposure Point 1.75 1.82 1.891.95 2.00 2.06 2.11 2.16

Note from the table that the tracking minimum for a single N54 magnet isabout eight radii from the center of the magnet to the nearest sensor(seven radii if using a weaker N24 magnet). Thus, the tracking minimumcan be a significant constraint on a magnet tracking system, andconsideration can be given to the sensor geometry required to meet thedesign requirements of a given application.

The worst-case scenario given above only remains valid when tracking nomore than a single magnet. The tracking boundaries are further limitedwhen tracking multiple magnets. To get the worst-case minimum distancefor multiple magnets, B_(max) can be divided by the number of magnetsbeing tracked. Table 3.2 lists the minimum distances from a pair ofmagnets to a magnetic field sensor.

TABLE 3.2 Minimum Magnet Distance from a Magnet Pair to an LSM9DS1 orLIS3MDL Sensor in Multiples of Magnet Radius, Given Residual FluxDensity Residual Flux Density (T) 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5(N-Rating) (N12) (N18) (N24) (N30) (N36) (N42) (N48) (N54) TrackingMinimum 8.25 8.58 8.89 9.18 9.45 9.70 9.95 10.18 Recalibration Point5.98 6.22 6.44 6.65 6.84 7.03 7.21 7.37 Maximum Exposure Point 2.21 2.292.38 2.45 2.52 2.59 2.66 2.72

Note that the tracking minimum distance for two magnets is about tworadii larger than the tracking minimum for one magnet. Again, this canbe considered a worst-case scenario because it assumes that the twomagnets occupy precisely the same location, which is not physicallypossible for magnets of non-zero diameter (though it can be a validassumption in an application where two magnets are located on oppositesides of a sensor). When a minimum separation between magnets isenforced in an application, this can be taken into account in the systemdesign when specifying the tracking boundaries.

Finally, though magnetic field disturbances can be compensated foralgorithmically, they can still add to the total magnetic field seen bythe sensors. Expected maximum magnetic field disturbance levels can beaccounted for in these tracking boundary calculations. This can beperformed by subtracting the expected maximum magnetic field disturbancefrom B_(max) before dividing it by the total number of magnets beingtracked. If full magnetic shielding is used, this may not be an issue.

The tables above rely on accurate radius and residual flux densityinformation. Deviations in magnet radius or magnet residual flux densitydue to manufacturing tolerances can necessitate additional bounding ofthe tracking volume.

The values in Tables 3.1 and 3.2 depend on which sensor is used. Forexample, using the LIS2MDL, the tracking minimum corresponds to a largerfull-scale range of 4.9152 mT, which is approximately the same as therecalibration point of the LIS3MDL. Furthermore, the LIS2MDL has amaximum exposure point corresponding to 1 Tesla, a magnetic field limitten times larger than the LIS3MDL can withstand (the magneticdisturbance field is not listed for the LIS2MDL, so the recalibrationpoint cannot be determined from the datasheet). The LIS2MDL, however,has a much lower sampling rate and suffers from increased nonlinearity,and these are also among the sensor characteristics that should beconsidered in the sensor selection for a particular application.

4. Detailed Explanation of Calibration

The following section presents improvements to the methods and systemsrelating to tracking of magnetic objects found in WO2019/074950, “Methodfor Neuromechanical and Neuroelectromagnetic Mitigation of LimbPathology;” the teachings of which are incorporated by reference hereinin their entirety.

A method of calibrating a magnetic field sensor array that comprises aplurality of magnets includes three-dimensionally rotating a magneticfield sensor array in a uniform magnetic field and recording data fromeach of the plurality of sensors. The method further includescalculating a non-rotating transformation of the recorded data using anellipsoid fit of each of the plurality of sensors and calculating ascaling factor of each of the plurality of sensors based on relativedimensions of the ellipsoid fit. A transformed, scaled dataset for eachof the plurality of sensors is obtained through application of thecalculated non-rotating transformation and calculated scaling factor.The obtained transformed, scaled datasets are rotated to align, therebyproviding for a determination of a relative orientation of each of theplurality of sensors to calibrate the magnetic field sensor array.

Magnetic field sensors are calibrated to account for “hard” and “soft”iron effects. “Hard” iron effects are offsets caused by magnetizedcomponents affixed to the same geometric frame as the sensor, and “soft”iron effects are magnetic field redirections caused by temporarilymagnetizable components affixed to the same geometric frame as thesensor. When multiple sensors are calibrated simultaneously, the scalingof the sensors relative to one another is also considered, as well astheir relative positions and orientations. In this section, methods ofaddressing the hard and soft iron effects are first provided. Then,recognizing the importance of position and orientation calibration theframework is extended to address how to perform position and orientationcalibration on one or more arrays of magnetic field sensors. Previouswork has relied upon fabrication tolerances for position and orientationalignment.

Throughout this section, the use of three-axis magnetic field sensors isassumed. Modifications to these methods can be made for single- ortwo-axis sensors.

4.1 Distortion Correction Via Ellipsoid Fitting with Relative Gains

Hard and soft iron effects can be corrected for on a magnetic fieldsensor by first recording the magnetic field data at the sensor whilerotating the sensor in a spatially-uniform ambient magnetic field (seeFIG. 22 ), and then transforming this ellipsoidally-distributed magneticfield data into a sphere centered at the origin. Though it is possibleto use just the midpoint between the minimum and maximum magnetic fielddata points along each axis of the rotated-sensor data to compensate forhard iron offsets, the ellipsoid-to-sphere transformation allowscompensation for soft iron effects. An ellipsoid-to-spheretransformation is also less prone to error when the rotated-sensordataset does not fully cover all possible sensor array orientations.Further, the ellipsoid-to-sphere calibration method has the advantage ofaccounting for any intrasensor coordinate axis non-orthogonalities dueto manufacturing misalignment.

In this section, ellipsoid-to-sphere calibration methods are provided. Astraightforward method for scaling the sensor gains by transforming therotated-sensor-array data from all sensors into a sphere of equivalentsize for the calibration of multiple sensors is also provided. With thetransformation parameters known, these parameters can then be applied toall future measured magnetic field data, even when the measured field isno longer spatially uniform (e.g., when tracking magnets).

An overview of a calibration method is shown in FIG. 24 . The method 900includes rotating a sensor array in three dimensions in a uniformmagnetic field (902) and calculating a non-rotating transformation usingan ellipsoid fit of the data set recorded by each sensor (904). Ascaling of the sensors from the relative dimensions of the ellipsoidfits is calculated (906). The transformed, scaled datasets are rotateduntil they align optimally with one another (908), and, using theparameters from the ellipsoid fits, scaling, rotation, sensor data ismodified while tracking is performed (910).

To ensure that the field is spatially uniform in the volume used forcalibration, calibration can be performed at a distance away frommagnets, walls, tables, chairs, floor, etc., such that these items donot affect the magnetic field sampled during the calibration.

The mathematics for computing and implementing the ellipsoid-to-spheretransformation for a magnetic field sensor array are provided below.Following the work from [20], for the ith sensor we seek a symmetricpositive-definite matrix

$\begin{matrix}{{A_{1} = \begin{bmatrix}A_{i} & {D_{i}/2} & {E_{i}/2} \\{D_{i}/2} & B_{i} & {F_{i}/2} \\{E_{i}/2} & {F_{i}/2} & C_{i}\end{bmatrix}},} & (4.1)\end{matrix}$

a vector b_(i)=[G_(i),H_(i),K_(i)]^(T), and a constant L_(i) such that,for {tilde over (B)}_(ik)=[{umlaut over (B)}_(ikx),{umlaut over(B)}_(iky),{umlaut over (B)}_(ikz)]^(T).

B _(ik) ^(T) A _(j) B _(ik) +b _(j) ^(T) B _(ik) +L _(i)=0  (4.2)

for all {B _(ik)}_(k=1) ^(K) in the calibration data set. This is theparametric equation of an arbitrary ellipsoid and can be expanded as:

$\begin{matrix}{{{\frac{1}{3}\left( {{\left( {A_{i} + B_{i} + C_{i}} \right)\left( {{\overset{\_}{B}}_{ikx}^{2} + {\overset{\_}{B}}_{iky}^{2} + {\overset{\_}{B}}_{ikz}^{2}} \right)} + {\left( {A_{i} - C_{i}} \right)\left( {{\overset{\_}{B}}_{ikx}^{2} + {\overset{\_}{B}}_{iky}^{2} + {\overset{\_}{B}}_{ikz}^{2}} \right)} + {\left( {A_{i} - B_{i}} \right)\left( {{\overset{\_}{B}}_{ikx}^{2} - {2{\overset{\_}{B}}_{iky}^{2}} + {\overset{\_}{B}}_{ikz}^{2}} \right)}} \right)} + {D_{i}{\overset{\_}{B}}_{ikx}{\overset{\_}{B}}_{iky}} + {E_{i}{\overset{\_}{B}}_{ikx}{\overset{\_}{B}}_{ikz}} + {F_{i}{\overset{\_}{B}}_{iky}{\overset{\_}{B}}_{ikz}} + {G_{i}{\overset{\_}{B}}_{ikx}} + {H_{i}{\overset{\_}{B}}_{iky}} + {K_{i}{\overset{\_}{B}}_{ikz}} + L_{i}} = 0} & (4.3)\end{matrix}$

Performing the following linear combinations,

A _(i) +B _(i) +C _(i) =W _(i)  (4.4a)

A _(i)−C_(i) =−U _(i) W _(i)  (4.4b)

A _(i) −B _(i) =−V _(i) W _(i)  (4.4c)

D _(i)=−4M _(i) W _(i)/3  (4.4d)

E _(i)=−2N _(i) W _(i)/3  (4.4e)

F _(i)=−2p _(i) W _(i)/3  (4.4f)

G _(i) =−Q _(i) W _(i)/3  (4.4g)

H _(i) =−R _(i) W _(i)/3  (4.4h)

K _(i) =−S _(i) W _(i)/3  (4.4i)

L _(i) =−T _(i) W _(i)/3  (4.4j)

allows (4.3) to be written as

$\begin{matrix}{{\frac{1}{3}\left( {{W_{i}\left( {{\overset{\_}{B}}_{ikx}^{2} + {\overset{\_}{B}}_{iky}^{2} + {\overset{\_}{B}}_{ikz}^{2}} \right)} - {U_{i}{W_{i}\left( {{\overset{\_}{B}}_{ikx}^{2} + {\overset{\_}{B}}_{iky}^{2} - {2{\overset{\_}{B}}_{ikz}^{2}}} \right)}} - {V_{i}{W_{i}\left( {{\overset{\_}{B}}_{ikx}^{2} - {2{\overset{\_}{B}}_{iky}^{2}} + {\overset{\_}{B}}_{ikz}^{2}} \right)}} - {4M_{i}W_{i}{\overset{\_}{B}}_{ikx}{\overset{\_}{B}}_{iky}} - {2N_{i}W_{i}{\overset{\_}{B}}_{ikx}{\overset{\_}{B}}_{ikz}} - {2P_{i}W_{i}{\overset{\_}{B}}_{iky}{\overset{\_}{B}}_{ikz}} - {Q_{i}W_{i}{\overset{\_}{B}}_{ikx}} - {R_{i}W_{i}{\overset{\_}{B}}_{iky}} - {S_{i}W_{i}{\overset{\_}{B}}_{ikz}} - {T_{i}W_{i}}} \right)} = 0} & (4.5)\end{matrix}$

Multiplying both sides by −3/W, and moving the leftmost term to theright-hand side then gives:

$\begin{matrix}{{{\left( {{\overset{\_}{B}}_{ikx}^{2} + {\overset{\_}{B}}_{iky}^{2} - {2{\overset{\_}{B}}_{ikz}^{2}}} \right)U_{i}} + {\left( {{\overset{\_}{B}}_{ikx}^{2} - {2{\overset{\_}{B}}_{iky}^{2}} + {\overset{\_}{B}}_{ikz}^{2}} \right)V_{i}} + {4{\overset{\_}{B}}_{ikx}{\overset{\_}{B}}_{iky}M_{i}} + {2{\overset{\_}{B}}_{ikx}{\overset{\_}{B}}_{ikz}N_{i}} + {2{\overset{\_}{B}}_{iky}{\overset{\_}{B}}_{ikz}P_{i}} + {{\overset{\_}{B}}_{ikx}Q_{i}} + {{\overset{\_}{B}}_{iky}R_{i}} + {{\overset{\_}{B}}_{ikz}S_{i}} + T_{i}} = {\left( {{\overset{\_}{B}}_{ikx}^{2} + {\overset{\_}{B}}_{iky}^{2} + {\overset{\_}{B}}_{ikz}^{2}} \right).}} & (4.6)\end{matrix}$

The variables in this equation can then be determined by usingleast-squares to solve

the matrix equation:

$\begin{matrix}{\text{?}} & (4.7)\end{matrix}$ ?indicates text missing or illegible when filed

or, most succinctly,

A_(i)s_(LS) _(i) =e_(i)  (4.8)

Once this equation is solved for the values of s_(LSi), these variablesmust be mapped back through the inverse of the linear combinations(4.4a)-(4.4j). First, we solve for A_(i), B_(i), and C_(i). From(6.4a)-(6.4c), we have

$\begin{matrix}{{\begin{bmatrix}1 & 1 & 1 \\1 & 0 & {- 1} \\1 & {- 1} & 0\end{bmatrix}\begin{bmatrix}A_{i} \\B_{i} \\C_{i}\end{bmatrix}} = \begin{bmatrix}W_{i} \\{{- U_{i}}W_{i}} \\{{- V_{i}}W_{i}}\end{bmatrix}} & (4.9)\end{matrix}$

so A_(i), B_(i), and C_(i) are given by

$\begin{matrix}\begin{matrix}{\begin{bmatrix}A_{i} \\B_{i} \\C_{i}\end{bmatrix} = {\begin{bmatrix}1 & 1 & 1 \\1 & 0 & {- 1} \\1 & {- 1} & 0\end{bmatrix}^{- 1}\begin{bmatrix}W_{i} \\{{- U_{i}}W_{i}} \\{{- V_{i}}W_{i}}\end{bmatrix}}} \\{= {{\frac{1}{3}\begin{bmatrix}1 & 1 & 1 \\1 & 1 & {- 2} \\1 & {- 2} & 1\end{bmatrix}}\begin{bmatrix}W_{i} \\{{- U_{i}}W_{i}} \\{{- V_{i}}W_{i}}\end{bmatrix}}} \\{= {{\frac{W_{i}}{3}\begin{bmatrix}1 & {- 1} & {- 1} \\1 & {- 1} & 2 \\1 & 2 & {- 1}\end{bmatrix}}\begin{bmatrix}1 \\U_{i} \\V_{i}\end{bmatrix}}}\end{matrix} & (4.1)\end{matrix}$

Now as described in [13], because we are transforming the data to asphere, we wish to have tr(A_(i))=3, or A_(i)+B_(i)+C_(i)=3, so we haveW_(i)=3, giving

$\begin{matrix}{\begin{bmatrix}A_{i} \\B_{i} \\C_{i}\end{bmatrix} = {\begin{bmatrix}1 & {- 1} & {- 1} \\1 & {- 1} & 2 \\1 & 2 & {- 1}\end{bmatrix}\begin{bmatrix}1 \\U_{i} \\V_{i}\end{bmatrix}}} & (4.11)\end{matrix}$

From (4.4d)-(4.4f), Di, Ei, and Fi are then given simply by

D _(i)=−4M _(i)  (4.12a)

E _(i)=−2N _(i)  (4.12b)

F _(i)=−2p _(i),  (4.12c)

and substituting (4.11)-(4.12c) into (4.1) gives

$\begin{matrix}{A_{i} = \begin{bmatrix}{1 - U_{i} - V_{i}} & {{- 2}M_{i}} & {- N_{i}} \\{{- 2}M_{i}} & {1 - U_{i} + {2V_{i}}} & {- P_{i}} \\{- N_{i}} & {- P_{i}} & {1 + {2U_{i}} - V_{i}}\end{bmatrix}} & (4.13)\end{matrix}$

Finally, we note from (4.4j) that the constant L_(i)=−T_(i).

To use the results of the least-squares fit (4.8) in transforming futuremagnetic field measurements, we first convert from the ellipsoidrepresentation of (4.2), the parametric form

B _(ik) ^(T) A _(i) B _(ik) +b _(i) ^(T) B _(ik) +L _(i)=0,  (4.2revisited)

to an ellipsoid representation where an offset is applied first,followed by a transformation, the quadric form

( B _(jk) −v _(i))^(T) A _(i)( B _(ik) −v _(i))=c _(i),  (4.14)

To represent the ellipsoid in this form, we first modify (4.14) so thatwe can calculate v_(i) and c_(i) by comparison with (4.2). Because A_(i)is symmetric (see (4.1)), we can expand (4.14) as

B _(ik) ^(T) A _(i) B _(ik) −B _(ik) ^(T) A _(i) v _(i) −v _(i) −v _(i)^(T) A _(i) B _(ik) +v _(i) ^(T) A _(i) v _(i) −c _(i)=0   (4.15a )

B _(ik) ^(T) A _(i) B _(ik)−( B _(ik) ^(T) A _(i) v _(i))^(T) −v _(i)^(T) A _(i) B _(ik) +v _(i) ^(T) A _(i) v _(i) −c _(i)=0   (4.15b)

B _(ik) ^(T) A _(i) B _(ik) −v _(i) ^(T) A _(i) ^(T) B _(ik) −v _(i)^(T) A _(i) B _(ik) +v _(i) ^(T) A _(i) v _(i) −c _(i)=0  (4.15c)

B _(ik) ^(T) A _(i) B _(ik)−2v _(i) ^(T) A _(i) B _(ik) +v _(i) ^(T) A_(i) v _(i) −c _(i)=0  (4.15d)

With (4.15d) and (4.2) in similar form, we can now equate coefficientsof like terms to get

$\begin{matrix}{{{- 2}v_{i}^{T}A_{i}} = b_{i}^{T}} & \left( {4.16a} \right)\end{matrix}$ $\begin{matrix}{\left( {{- 2}v_{i}^{T}A_{i}} \right)^{T} = \left( b_{i}^{T} \right)^{T}} & \left( {4.16b} \right)\end{matrix}$ $\begin{matrix}{{{- 2}A_{i}v_{i}} = b_{i}} & \left( {4.16c} \right)\end{matrix}$ $\begin{matrix}{v_{i} = {{- \frac{1}{2}}A_{i}^{- 1}b_{i}}} & \left( {4.16d} \right)\end{matrix}$ $\begin{matrix}{v_{i} = {{- \frac{1}{2}}{A_{i}^{- 1}\begin{bmatrix}G_{i} \\H_{i} \\K_{i}\end{bmatrix}}}} & \left( {4.16e} \right)\end{matrix}$ $\begin{matrix}{v_{i} = {{- \frac{1}{2}}{A_{i}^{- 1}\begin{bmatrix}Q_{i} \\R_{i} \\S_{i}\end{bmatrix}}}} & \left( {4.16f} \right)\end{matrix}$

where the vector vi contains the coordinates of the center of theellipsoid, which are used to compensate for hard iron effects byoffsetting the data. Equating the constant part of (4.14) with theconstant part of (4.6) gives

L _(i) =V _(i) ^(T) A _(i) v _(i) −c _(i)  (4.17a)

c _(i) =T _(i) +v _(i) ^(T) A _(i) v _(i)  (4.17b)

which we will use to scale the sensor measurements relative to oneanother. From Eqns. 11 and 12 of [11], this constant c_(i) isrepresentative of the square of the field strength seen by the sensor(i.e., the constant c_(i) is equal to the square of the radius of thesphere after sensitivity scaling has been performed, or c_(i)=r_(i) ²).We thus take the square root of this constant to get the sphere radius,r_(i)=√c_(i). Using the average of the sphere radii across all Nsensors,

$\begin{matrix}{{r_{avg} = {\frac{1}{N}{\sum\limits_{i = 1}^{N}r_{i}}}},} & (4.18)\end{matrix}$

a scaling factor can then be computed for each of the sensors,

$\begin{matrix}{g_{i} = {\frac{r_{avg}}{r_{i}}.}} & (4.19)\end{matrix}$

Using the method of [11] to perform a non-rotating transformation (sothat the sensors do not rotate relative to one another), we compute theinverse soft-iron matrix as the square root of the ellipsoid matrix

W _(i) ⁻¹ =A _(i) ^(½).  (4.20)

Once the calibration is performed, then while collecting data, for alldata streams coming in, the data is offset, transformed, and scaled viathe equation

{hacek over (B)} _(i) =g _(i) S _(i) W _(i) ⁻¹({tilde over (B)} _(i) −v_(i)),  (4.21)

where S_(i) is the sensor sensitivity given in the sensor's datasheet.

Performing this ellipsoid-to-sphere calibration and applying (4.21) tothe calibration data gives a set of spheres which are approximatelyequal and concentric, as shown in FIGS. 23A-C.

4.2 Extension to Arbitrarily Oriented Magnetic Field Sensors

If the magnetic field sensors are not aligned along the same axes (e.g.,if they are affixed to a curved rigid surface or to multiple,rigidly-connected but unaligned flat surfaces), the sensor orientationsrelative to one another (and, importantly, relative to the globalmagnetic field sensor array coordinate system) can be determined.

To calibrate these relative orientations, the point clouds from thedistortion correction calibration, after offset, transformation, andscaling, can be rotated until they are optimally aligned with oneanother. Because the calibration data is ordered (points are matched ateach time point), this alignment can be performed on a point-wise basisinstead of resorting to more complex point set registration techniques.

We first define the measured magnetic field B^({tilde over ( )}) _(i) asa rotated version of the sensor data after the distortion correctioncalibration.

$\begin{matrix} & (4.22)\end{matrix}$ $\begin{matrix}{{\overset{\sim}{B}}_{i}\overset{\Delta}{=}{{R_{z}\left( \alpha_{i} \right)}{R_{y}\left( \beta_{i} \right)}{R_{x}\left( \gamma_{i} \right)}{\overset{˜}{B}}_{i}}} \\{= {{{\begin{bmatrix}{\cos\alpha_{i}} & {{- s}{in}\alpha_{i}} & 0 \\{\sin\alpha_{i}} & {\cos\alpha_{i}} & 0 \\0 & 0 & 1\end{bmatrix}\begin{bmatrix}{\cos\beta_{i}} & 0 & {\sin\beta_{i}} \\0 & 1 & 0 \\{{- \sin}\beta_{i}} & 0 & {\cos\beta_{i}}\end{bmatrix}}\begin{bmatrix}1 & 0 & 0 \\0 & {\cos\gamma_{i}} & {{- s}{in}\gamma_{i}} \\0 & {\sin\gamma_{i}} & {\cos\gamma_{i}}\end{bmatrix}}\begin{bmatrix}{\check{B}}_{ix} \\{\check{B}}_{iy} \\{\check{B}}_{iz}\end{bmatrix}}} \\{= {{\begin{bmatrix}{\cos\alpha_{i}} & {{- s}{in}\alpha_{i}} & 0 \\{\sin\alpha_{i}} & {\cos\alpha_{i}} & 0 \\0 & 0 & 1\end{bmatrix}\begin{bmatrix}{\cos\beta_{i}} & {\sin\beta_{i}\sin\gamma_{i}} & {\sin\beta_{i}\cos\gamma_{i}} \\0 & {\cos\gamma_{i}} & {{- s}{in}\gamma_{i}} \\{{- s}{in}\beta_{i}} & {\cos\beta_{i}\sin\gamma_{i}} & {\cos\beta_{i}\cos\gamma_{i}}\end{bmatrix}}\begin{bmatrix}{\check{B}}_{ix} \\{\check{B}}_{iy} \\{\check{B}}_{iz}\end{bmatrix}}} \\{= {\begin{bmatrix}{\cos\alpha_{i}\cos\beta_{i}} & \begin{matrix}{{\cos\alpha_{i}\sin\beta_{i}\sin\gamma_{i}} -} \\{\sin\alpha_{i}\cos\gamma_{i}}\end{matrix} & \begin{matrix}{{\cos\alpha_{i}\sin\beta_{i}\cos\gamma_{i}} +} \\{\sin\alpha_{i}\sin\gamma_{i}}\end{matrix} \\{\sin\alpha_{i}\cos\beta_{i}} & \begin{matrix}{{\sin\alpha_{i}\sin\beta_{i}\sin\gamma_{i}} +} \\{\cos\alpha_{i}\cos\gamma_{i}}\end{matrix} & \begin{matrix}{{\sin\alpha_{i}\sin\beta_{i}\cos\gamma_{i}} -} \\{\cos\alpha_{i}\sin\gamma_{i}}\end{matrix} \\{{- \sin}\beta_{i}} & {\cos\beta_{i}\sin\gamma_{i}} & {\cos\beta_{i}\cos\gamma_{i}}\end{bmatrix}\begin{bmatrix}{\check{B}}_{ix} \\{\check{B}}_{iy} \\{\check{B}}_{iz}\end{bmatrix}}}\end{matrix},$

where α_(i), β_(i), γ_(i) are respectively our estimates of the yaw,pitch, and roll of the ith sensor. When no rotation is required (e.g.,if the sensors are precisely aligned during manufacturing and assembly),B^({tilde over ( )}) _(i)=B^({hacek over ( )}) _(i). A reference sensoris then chosen, which we will denote as sensor n. This reference sensorwill define the orientation of the axes of the sensor array coordinatesystem. This sensor is defined to be oriented, such thatα_(n)=β_(n)=γ_(n)=0 (and, as a consequence, B^({tilde over ( )})_(n)=B^({hacek over ( )}) _(n)). The rotation parameters of theremaining sensors from the reference sensor are then each given by thesolution to the optimization

$\begin{matrix}{{\underset{\alpha_{i},\beta_{i},\gamma_{i}}{argmin}{\sum\limits_{k = 1}^{K}{{{\overset{\sim}{B}}_{nk} - {R_{i}{\check{B}}_{ik}}}}_{2}}},} & (4.23)\end{matrix}$

where we have defined R_(i), R_(z)(α_(i))R_(y)(β_(i))R_(x)(γ_(i)). Uponclose inspection of (4.23), it is found that the structure of thisproblem has previously received particular attention in the Kabschalgorithm (as a special case of the orthogonal Procrustes problem, inturn a special case of Wahba's problem) and can be solved using singularvalue decomposition [4, 5].

Once the rotation matrix R_(i) for each sensor is computed, the measuredmagnetic field value (the final value used for comparison with predictedmagnetic field values in the optimization) then becomes

{tilde over (B)} _(i) =R _(i) g _(i) S _(i) W _(i) ⁻¹( B _(i) −v_(i)).  (4.24)

Defining T_(i)=R_(i)g_(i)S_(i)W_(i) ⁻¹ gives

{tilde over (B)} _(i) =T _(i)( B _(i) −v _(i)),  (4.25)

or augmenting T_(i) with −v_(i) to get the affine transformation matrixT_(i)′=[R_(i)g_(i)SW_(i) ⁻¹|−v_(i)], and using the augmented 4-vector B_(ik)′=[B _(ikx), B _(iky), B _(ikz), 1]^(T), we have simply

{tilde over (B)} _(i) =T _(i) ′B _(i)′.  (4.26)

4.3 Extension to Arbitrarily Positioned Magnetic Field Sensors

In Section 3 herein, we assumed that the positions of the magnetic fieldsensors were previously known while measuring the dipole strength of oneor more magnets. Below, we present a method for determining sensorpositions using a spherical magnet for position calibration, extendingthe technique of Section 3.3.

If the sensor positions are unknown, this calibration can simultaneouslycompute the sensor positions and the dipole's strength. This can benecessary because a magnet's strength is not consistent as viewed by onesensor array versus another.

As before, the magnetic field prediction error is used for theoptimization cost function. At the ith sensor and kth timestep, themagnetic field prediction error, E_(ik), is the difference between themeasured magnetic field B^({tilde over ( )}) _(ik) and the predictedmagnetic field B_(ik),

E _(ik) =B _(ik) −{tilde over (B)} _(ik).  (4.27)

Equation (4.27) is identical to (3.4), with the difference being in theway that B_(ik) is calculated, as a function of sensor positions.

For each timestep (the kth timestep), the derivatives of the magneticfield prediction errors with respect to each of the sensor positionestimates, s_(ix), s_(iy), and s_(iz), are calculated and placed in aseparate submatrix. Because the measured magnetic fieldB^({tilde over ( )}) _(ik) of measurement k does not vary with respectto the estimated sensor positions (for example,∂/∂s_(xi)E_(ik)=∂/∂s_(ix)B_(ik)), the derivatives of the cost functioncan be written as a function of B_(ik) alone. Thus, the Jacobiansubmatrix for sensor position error corresponding to the ith sensor andkth measurement can be calculated as

$\begin{matrix}\begin{matrix}{P_{ik} = \begin{bmatrix}{\frac{\partial}{\partial s_{ix}}E_{ikx}} & {\frac{\partial}{\partial s_{iy}}E_{ikx}} & {\frac{\partial}{\partial s_{iz}}E_{ikx}} \\{\frac{\partial}{\partial s_{ix}}E_{iky}} & {\frac{\partial}{\partial s_{iy}}E_{iky}} & {\frac{\partial}{\partial s_{iz}}E_{iky}} \\{\frac{\partial}{\partial s_{ix}}E_{ikz}} & {\frac{\partial}{\partial s_{iy}}E_{ikz}} & {\frac{\partial}{\partial s_{iz}}E_{ikz}}\end{bmatrix}} \\{= \begin{bmatrix}{\frac{\partial}{\partial s_{ix}}B_{ikx}} & {\frac{\partial}{\partial s_{iy}}B_{ikx}} & {\frac{\partial}{\partial s_{iz}}B_{ikx}} \\{\frac{\partial}{\partial s_{ix}}B_{iky}} & {\frac{\partial}{\partial s_{iy}}B_{iky}} & {\frac{\partial}{\partial s_{iz}}B_{iky}} \\{\frac{\partial}{\partial s_{ix}}B_{ikz}} & {\frac{\partial}{\partial s_{iy}}B_{ikz}} & {\frac{\partial}{\partial s_{iz}}B_{ikz}}\end{bmatrix}}\end{matrix} & (4.28)\end{matrix}$

The columns of derivatives in (4.28) are written in terms of s_(ix),s_(iy), and s_(iz), but B_(i) in Eqns. 8 a and 8 c of the TaylorReference is written in terms of x⁻ _(ij), y⁻ _(ij), and z⁻ _(ij). Fromthe definitions of x⁻ _(ij), y⁻ _(ij), and z⁻ _(ij) in (2.3) we see thatx⁻ _(ij), y⁻ _(ij), and z⁻ _(ij) are functions of s_(ix), s_(iy), ands_(iz), respectively. The chain rule can thus be used for all M magnetsat the kth measurement as

$\begin{matrix}{\frac{\partial B_{ik}}{\partial s_{ix}} = {{\sum\limits_{j^{\prime} = 1}^{M}\left( {{\frac{\partial B_{ik}}{\partial{\overset{\_}{x}}_{{ij}^{\prime}k}}\frac{\partial{\overset{\_}{x}}_{{ij}^{\prime}k}}{\partial s_{ix}}} + {\frac{\partial B_{ik}}{\partial{\overset{\_}{y}}_{{ij}^{\prime}k}}\frac{\partial{\overset{\_}{y}}_{{ij}^{\prime}k}}{\partial s_{ix}}} + {\frac{\partial B_{ik}}{\partial{\overset{\_}{z}}_{{ij}^{\prime}k}}\frac{\partial{\overset{\_}{z}}_{{ij}^{\prime}k}}{\partial s_{ix}}}} \right)} = {\sum\limits_{j^{\prime} = 1}^{M}\left( {{\frac{\partial B_{ik}}{\partial{\overset{\_}{x}}_{{ij}^{\prime}k}}(1)} + {\frac{\partial B_{ik}}{\partial{\overset{\_}{y}}_{{ij}^{\prime}k}}(0)} + {\frac{\partial B_{ik}}{\partial{\overset{\_}{z}}_{{ij}^{\prime}k}}(0)}} \right)}}} & (4.29)\end{matrix}$

for s_(ix), and similarly for s_(iy) and s_(iz) to get

$\begin{matrix}{P_{ij} = {\sum\limits_{j^{\prime} = 1}^{M}\text{?}}} & (4.3)\end{matrix}$ ?indicates text missing or illegible when filed

in which we note that the derivative for each element is now computed asa summation over all M terms of Eqn. 8 (of the Taylor Reference). Inother words, each of the elements of P_(ik) is a negated summation ofthe M elements of J corresponding to the ith sensor and kth measurementalong each given axis. However, unless some advantage is found forposition calibration using multiple magnets, the use of a single magnet(M=1) can be preferred, and in this case no summation is required.

The submatrix of the Jacobian of magnetic field prediction errorscorresponding to all sensor position estimates at measurement k is thengiven by

$\begin{matrix}{P_{k} = {\begin{bmatrix}P_{1k} & 0 & \ldots & 0 \\0 & P_{2k} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots \\0 & 0 & \ldots & P_{Nk}\end{bmatrix}.}} & (4.31)\end{matrix}$

The dipole measurement matrix can then be augmented to give us thesensor position calibration matrix

$\begin{matrix}{\overset{\sim}{\overset{\sim}{J}} = \begin{bmatrix}J_{1} & 0 & \ldots & 0 & M_{1} & P_{1} \\0 & J_{2} & \ldots & 0 & M_{2} & P_{2} \\ \vdots & \vdots & \ddots & & \vdots & \vdots \\0 & 0 & \ldots & J_{K} & M_{K} & P_{K}\end{bmatrix}} & (4.32)\end{matrix}$

where K is the number of samples taken at different time steps.

Note that this derivation above does not tell a perfect picture, sinceeither a few distinct positions must be known (or some combination ofknown position and/or orientation for multiple points), or one of thesensors can be used as the “origin” sensor (likely a simpler option).Say that the reference sensor is sensor n. Then, P_(k) is given by

$\begin{matrix}{P_{k} = \begin{bmatrix}P_{1k} & 0 & \ldots & 0 & 0 & \ldots & 0 \\0 & P_{2k} & \ldots & 0 & 0 & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\0 & 0 & \ldots & P_{{({n - 1})}k} & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & 0 & \ldots & 0 \\0 & 0 & \ldots & 0 & P_{{({n + 1})}k} & \ldots & 0 \\ \vdots & \vdots & \ddots & \vdots & \vdots & \ddots & \vdots \\0 & 0 & \ldots & 0 & 0 & \ldots & P_{Nk}\end{bmatrix}} & (4.33)\end{matrix}$

Using the notation given above and letting n=1, and when calibratingpositions with a single magnet, the sensor position calibration matrixappears in expanded form as shown in FIG. 36 .

4.4 Further Technical Considerations

Throughout this section, the assumption was made that the sensors arelinear. This assumption had no bearing on the given approach tocalibrating sensor orientations and positions, or on the given methodfor measuring dipole strengths. However, nonlinearity of the sensors(e.g., of the LIS3MDL or LSM9DS1 magnetic tunnel junctions) can warp themagnetic field measurements, and soft iron saturation can warp theactual magnetic field at the sensor. Though it was suggested in theTaylor Reference that the magnetization inhomogeneity of sphericalpermanent magnets can contribute to the tracking error observed whenmagnets are close to the sensing array, tracking errors at near-fieldcaused by sensor nonlinearity are likely greater in magnitude. Uponinspection of the datasheet for the LIS3MDL (which is equivalent to theLSM9DS1), when sensing a 600 μT magnetic field with a full-scale rangeof 1200 μT, the typical specification of the nonlinearity is 0.12% offull-scale range, or 1.44 μT, roughly half the RMS noise of the sensor.This error can be even more substantial when sensing larger magneticfields, as well as when sensing with a full-scale range of 1600 μT. Itmay thus be valuable to consider implementing nonlinearity calibrationon the sensors to increase tracking accuracy at near-field. However,such calibration techniques can be non-trivial to implement, and incombining and extending the strategies above, can potentially adverselyaffect tracking latency.

In the Taylor Reference, spike noise was filtered out using athree-point median filter. It was originally thought that the spikenoise was quantum mechanical shot noise from the magnetic tunneljunction of the LSM9DS1 tunneling magnetoresistors used in the TaylorReference. With further investigation, however, it was discovered thatthis spike noise was just due to communication issues, which wascorrected by enabling the block data update bit on the LIS3MDL tunnelingmagnetoresistors described in Example 1 herein.

Upon further inspection of the LIS3MDL datasheet, the LIS3MDL has atemperature sensor output change vs. temperature of 8 LSB/° C. Thus, asingle degree change in temperature can cause a magnetic field sensingerror of 464 nT, greater even than the RMS noise of 410 nT reported forthe sensor. To account for these temperature effects, as much aspossible the sensors were calibrated at the same temperature as wasexpected in all experiments. However, a swing in temperature of, forinstance, the ten degrees Celsius difference between room temperatureand the temperature of a turkey leg, can easily skew the magnetic fieldreadings by 11 times the RMS noise. If this skewing is normallydistributed, when using many sensors, this should have a negligibleeffect on the tracking accuracy in our application.

In the Taylor Reference a backwards labeling of the magnetic fieldsensing axes in the LSM9DS1 datasheet was erroneously followed. Because,mathematically, the reversal of all three axes is equivalent to thetracking of magnets that are reversed in polarity, this went undetecteduntil the seventh-generation sensor array employed in Example 1, when abar magnet was used to aid in setting up the sensor array coordinateaxes. Thus, in the Taylor Reference, the magnets were tracked as if theywere aligned south-pole-up in the tracking wands, instead ofnorth-pole-up as intended and physically implemented, much to ourconfusion during tracking, but without any effect on the validity of theresults.

5. Extension of Disturbance Rejection to Higher Order Terms of theDisturbance

As discussed in the Taylor Reference, magnetic disturbance fields can becompensated for by tracking the disturbance fields as part of a magnettracking optimization algorithm. When objects causing or warpingdisturbance fields are sufficiently far away, the disturbance presentsas a spatially-uniform field. For example, when moving a magneticcompass side-to-side far away from any magnets, the reading on thecompass changes very little, indicating that the field is relativelyconsistent within the local volume traversed. For example, FIG. 25 showsmagnetic field fluctuations as measured inside a train when being passedby another train traveling in the opposite direction on an adjacenttrack. Notice that the corresponding x, y, and z axes across all of thesensors move up and down together. In other words, a uniform disturbancefield assumption can reject all far-field magnetic disturbances,regardless of the source.

However, as these objects come nearer to the magnetic field sensorarray, the disturbance becomes increasingly non-uniform. That is, thedisturbance presents as a spatial gradient. Even in FIG. 25 , morecareful inspection reveals that the measurements across equivalent axesare not precisely equal. These nonuniform disturbances can becompensated for in much the same way that non-uniform disturbances werecompensated for in the Taylor Reference, by including the spatialgradient components of the disturbance field in the tracking parameters.

5.1 Compensation for the Spatial Gradient of the Disturbance

Modifying Eqn. 8 of the Taylor Reference to include the spatial gradientof the disturbance field gives

$\begin{matrix}{{B_{i} = {G + {s_{i} \cdot \mathcal{G}} + {\sum\limits_{j^{\prime} = 0}^{j^{\prime} = M}\left( {\frac{3{r_{{ij}^{\prime}}\left( {{\overset{\_}{m}}_{j^{\prime T}} \cdot r_{{ij}^{\prime}}} \right)}}{r_{{ij}^{\prime}}^{5}} - \frac{{\overset{\_}{m}}_{j^{\prime}}}{r_{{ij}^{\prime}}^{3}}} \right)}}},} & (5.1)\end{matrix}$

where G is the spatially-uniform (zeroth order) component of thedisturbance field and G is a tensor containing the first-order spatialgradient information. In expanded form, this is shown as

$\begin{matrix}{{B_{i} = {G + {\begin{bmatrix}s_{ix} & s_{iy} & s_{iz}\end{bmatrix}\begin{bmatrix}G_{x_{x}} & G_{y_{x}} & G_{z_{x}} \\G_{x_{y}} & G_{y_{y}} & G_{z_{y}} \\G_{x_{z}} & G_{y_{z}} & G_{z_{z}}\end{bmatrix}} + {\sum\limits_{j^{\prime} = 0}^{j^{\prime} = M}\left( {\frac{3{r_{{ij}^{\prime}}\left( {{\overset{\_}{m}}_{j^{\prime T}} \cdot r_{{ij}^{\prime}}} \right)}}{r_{{ij}^{\prime}}^{5}} - \frac{{\overset{\_}{m}}_{j^{\prime}}}{r_{{ij}^{\prime}}^{3}}} \right)}}},} & (5.2)\end{matrix}$

where G_(xy) is the first-spatial-derivative with respect toy of thex-directed disturbance field. Note that, fully generalizing, G_(xy) andG_(yx) represent distinct properties. Now G introduces nine additionalparameters to the optimization, so a new submatrix is concatenated tothe Jacobian matrix

$\begin{matrix}\begin{matrix}{{Dxyz}_{i} = \begin{bmatrix}{\frac{\partial}{\partial G_{x_{x}}}B_{ix}} & {\frac{\partial}{\partial G_{y_{x}}}B_{ix}} & {\frac{\partial}{\partial G_{z_{x}}}B_{ix}} & {\frac{\partial}{\partial G_{x_{y}}}B_{ix}} & {\frac{\partial}{\partial G_{y_{y}}}B_{ix}} & {\frac{\partial}{\partial G_{z_{y}}}B_{ix}} & {\frac{\partial}{\partial G_{x_{z}}}B_{ix}} & {\frac{\partial}{\partial G_{y_{z}}}B_{ix}} & {\frac{\partial}{\partial G_{z_{z}}}B_{ix}} \\{\frac{\partial}{\partial G_{x_{x}}}B_{iy}} & {\frac{\partial}{\partial G_{y_{x}}}B_{iy}} & {\frac{\partial}{\partial G_{z_{x}}}B_{iy}} & {\frac{\partial}{\partial G_{x_{y}}}B_{iy}} & {\frac{\partial}{\partial G_{y_{y}}}B_{iy}} & {\frac{\partial}{\partial G_{z_{y}}}B_{iy}} & {\frac{\partial}{\partial G_{x_{z}}}B_{iy}} & {\frac{\partial}{\partial G_{y_{z}}}B_{iy}} & {\frac{\partial}{\partial G_{z_{z}}}B_{iy}} \\{\frac{\partial}{\partial G_{x_{x}}}B_{iz}} & {\frac{\partial}{\partial G_{y_{x}}}B_{iz}} & {\frac{\partial}{\partial G_{z_{x}}}B_{iz}} & {\frac{\partial}{\partial G_{x_{y}}}B_{iz}} & {\frac{\partial}{\partial G_{y_{y}}}B_{iz}} & {\frac{\partial}{\partial G_{z_{y}}}B_{iz}} & {\frac{\partial}{\partial G_{x_{z}}}B_{iz}} & {\frac{\partial}{\partial G_{y_{z}}}B_{iz}} & {\frac{\partial}{\partial G_{z_{z}}}B_{iz}}\end{bmatrix}} \\{= \begin{bmatrix}s_{ix} & 0 & 0 & s_{iy} & 0 & 0 & s_{iz} & 0 & 0 \\0 & s_{ix} & 0 & 0 & s_{iy} & 0 & 0 & s_{iz} & 0 \\0 & 0 & s_{ix} & 0 & 0 & s_{iy} & 0 & 0 & s_{iz}\end{bmatrix}} \\{= \begin{bmatrix}{s_{ix}I_{3}} & {s_{iy}I_{3}} & {s_{iz}I_{3}}\end{bmatrix}}\end{matrix} & (5.3)\end{matrix}$

The further augmented Jacobian matrix including the gradient of thedisturbance field, which we'll call J^({circumflex over ( )}){circumflexover ( )}, is then given by

$\begin{matrix}{\hat{\overset{\hat{}}{J}} = \begin{bmatrix}J_{11} & J_{12} & \ldots & J_{1M} & I_{3} & {s_{1x}I_{3}} & {s_{1y}I_{3}} & {s_{1z}I_{3}} \\J_{21} & J_{22} & \ldots & J_{2M} & I_{3} & {s_{2x}I_{3}} & {s_{2y}I_{3}} & {s_{2z}I_{3}} \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots & \vdots & \vdots \\J_{N1} & J_{N2} & \ldots & J_{NM} & I_{3} & {s_{Nx}I_{3}} & {s_{Ny}I_{3}} & {s_{Nz}I_{3}}\end{bmatrix}} & (5.4)\end{matrix}$

Note that some of these elements are not meaningful when working with aplanar array, so they may not be used in practice when the sensors lieon a single flat circuit board.

One could continue to higher-order dimensions, with a 27-element tensorfor the second-order, an 81-element tensor for the third-order, andfurther, but the added disturbance tracking sophistication can degradetracking performance, especially when using a smaller number of sensors.

At some level, when an intricately complicated disturbance field isexpected, alternative solutions can be considered. For instance,additional magnetic dipoles can be tracked (with the earth itself beingtracked as a large dipole far away). Alternatively, lightweightstrategies for implementing magnetic shielding in ways that eliminate orsimplify the magnetic disturbance can be employed.

5.1.1 Simplification of the Spatial Gradient Disturbance Compensation

On the other hand, if we allow the assumption that the source of thedisturbance(s) is dipolar in nature (e.g., that the magnetic field isnot being generated by electric current in, for instance, a long wire),we can apply the conditions of symmetry and zero-trace to G, as in [9],to simplify the disturbance tracking, as

$\begin{matrix}{{B_{i} = {G_{0} + {\begin{bmatrix}s_{ix} & s_{iy} & s_{iz}\end{bmatrix}\begin{bmatrix}G_{x_{x}} & G_{x_{y}} & G_{z_{x}} \\G_{x_{y}} & G_{y_{y}} & G_{z_{y}} \\G_{z_{x}} & G_{z_{y}} & {{- G_{x_{x}}} - G_{y_{y}}}\end{bmatrix}} + {\sum\limits_{j^{\prime} = 0}^{j^{\prime} = M}\left( {\frac{3{r_{{ij}^{\prime}}\left( {{\overset{\_}{m}}_{j^{\prime T}} \cdot r_{{ij}^{\prime}}} \right)}}{r_{{ij}^{\prime}}^{5}} - \frac{{\overset{\_}{m}}_{j^{\prime}}}{r_{{ij}^{\prime}}^{3}}} \right)}}},} & (5.5)\end{matrix}$

where G can now be described by five parameters. We then have thereduced-size disturbance-gradient Jacobian submatrix

? ?indicates text missing or illegible when filed

which further simplifies the Jacobian matrix to

$\begin{matrix}{\hat{J} = {\begin{bmatrix}J_{11} & J_{12} & \ldots & J_{1M} & I_{3} & {\overset{\_}{D}}_{{xyz}_{1}} \\J_{21} & J_{22} & \ldots & J_{2M} & I_{3} & {\overset{\_}{D}}_{{xyz}_{2}} \\ \vdots & \vdots & \ddots & \vdots & \vdots & \vdots \\J_{N1} & J_{N2} & \ldots & J_{NM} & I_{3} & {\overset{\_}{D}}_{{xyz}_{N}}\end{bmatrix}.}} & (5.7)\end{matrix}$

5.2 Flexible Magnetic Shielding

In magnetomicrometry, as discussed above, uniform, linear, and evennonlinear magnetic disturbance fields can be compensated for usingsoftware. However, if the disturbance field is linear or nonlinear, thiscomplicates the software, and becomes increasing more difficult if thefield needs to be described by higher order terms. Further, if themagnetic disturbance is large, there is the potential for damage to themagnetic field sensors.

As discussed above, one solution to the disturbance problem is tointroduce magnetic shielding. Past use of the shielding technique fortracking magnetic beads in the body have used the full-room (Faradaycage) shielding, or placement of a cylindrical shield, open on bothsides, around the part of the body wherein magnetic beads are to betracked.

The full-room shielding is limited because it does not allow fordisturbance compensation outside the room, so it isn't a portablesolution. The cylindrical solution is limited because it must be largein order to fit the limb, and generalizing it requires creating asolution that is unnecessarily large for wrapping around, for instance,extremities that are smaller.

Devices and methods for shielding which surpasses these limitations andallows for items of variable geometries to be optimally shielded areprovided. Two examples of wearable shielding assemblies 1000, and 1050are described below. In both of these examples, the shielding isexternal to the magnet-sensor system. In either of these examples, theshielding can be embedded in the socket as part of a lamination orplaced around the outside of the socket.

As illustrated in FIGS. 26 and 27 , a wearable shielding assembly 1000,1050 can include an array of sensors (e.g., sensors 106 a-e, or magneticfield sensors 208 a-b, sensors not visible in FIGS. 26 and 27 )configured to detect a state change of at least magnetic targetimplanted at a tissue (e.g., targets 101, 102, 201, 202, targets notvisible in FIGS. 26 and 27 ). The shielding assembly further includes awearable receptacle 1010 within which the array of sensors can bedisposed. A geometrically-reconfigurable material 1020 is disposed aboutthe wearable or is integral with the wearable. Examples ofgeometrically-reconfigurable material are provided below. Thegeometrically-reconfigurable material can be flexible and can,optionally, be rigidly affixed to the wearable receptacle afterapplication.

5.2.1 Tiled Magnetic Shielding

FIG. 26 illustrates an example shielding 1000 that includes a pluralityof ferromagnetic tiles 1030. In this example, multiple ferromagnetictiles can be connected together to create a full magnetic shield. Theconnections between the multiple tiles in this invention can be made byconnectors 1032, such as, for example snap fasteners, hooks, buttons,hook-and-loop fasteners (such as Velcro, or 3M Dual-Lock fasteners),glue, suture threads, knots, and elastic bands. These connections can beestablished prior to use or can be adjusted by the user. Theferromagnetic tiles in this invention can be thin or thick and may bemulti-material (multi-layered), including MuMetal® (Magnetic ShieldCorp., IL), steel, and other ferromagnetic materials. These tiles can bea repeated size and shape but may also vary in shape and size to providefor a customizable fit.

Unlike in lamellar armor, the individual tiles, called plates or scalesin lamellar armor, need not overlap. Thus, the lacing or other means ofconnection between the scales can be made of a flexible material, withgaps being allowed between tiles when the shielding is stretched.Alternatively, the tiles may be attached to, for example, the outside ofa socket, using a method of attachment such as stickers, glue, epoxy orVelcro disposed on a back side of each of the tiles, in place of or inaddition to the connections directly between the tiles.

5.2.2 Layered Metal Mesh for Magnetic Shielding

FIG. 27 illustrates another example of shielding 1050 that includes aferromagnetic mesh 1040. In this example, metal mesh (also known ashardware cloth or chicken wire, depending upon material and context) isapplied around the region where magnets targets are being tracked.Application can occur as a single wrap around the region, for instance,around a lower leg as illustrated, or can occur as multiple wrappings ofone material or of multiple different materials. For example, one ormore layers of MuMetal mesh can be applied as a base layer for theshielding, followed by steel chicken wire for additional outer layers.Full coverage of the metal mesh is not necessary for magnetic shieldingto occur, so a small number of wraps can be used to minimize the weightand volume of the shielding. The wrap can be incorporated into thecomposite structure of a prosthetic socket, or incorporated into acomposite structure of an orthotic interface between an exoskeleton andthe human body.

6. Common Subexpression Elimination, Practical Notes, and FurtherImprovements

The following section presents improvements to the methods and systemsrelating to tracking of magnetic objects found in WO2019/074950, “Methodfor Neuromechanical and Neuroelectromagnetic Mitigation of LimbPathology;” the teachings of which are incorporated by reference hereinin their entirety.

A method of tracking one or more permanent magnets includes detecting asignal from each of the one or more permanent magnets, calculating ananalytically-derived Hessian matrix with respect to the detectedsignals, and determining a state of each of the one or more permanentmagnets based on the calculated Hessian matrix.

The method can further include calculating a predicted magnetic fieldvalue for each of the one or more permanent magnets and calculating aJacobian matrix with respect to the detected signals. Determining thestate of each of the one or more permanent magnets can be further basedon the calculated predicted magnetic field values and the calculatedJacobian matrix. Calculation of the Hessian matrix, the Jacobian matrix,and the predicted magnetic field values can be performed simultaneouslyusing common subexpression elimination. The calculation of the Hessianmatrix can be parallelized.

6.1 Common Subexpression Elimination

There are many common subexpressions encountered when calculating themagnetic field prediction error and its derivatives. As discussed in theTaylor Reference, to increase the speed of the algorithm, commonsubexpression elimination was implemented, evaluating repeatedsubexpressions only once whenever possible. Noting that theimplementation of common subexpression elimination is besthand-optimized [15], the common subexpression elimination solution usedwas first developed on paper before being implemented on a computer.Local variables were used in the code whenever possible, and the numberof variables used was kept as low as possible to increase the use offast-access registers for memory storage (a technique most efficientlyimplemented in a compiled language). This section describes thestructure of subexpression calculations that can be used to achievelow-latency tracking.

In the equations below, some expressions represent variabledeclarations, some expressions indicate self-referential variablemodification, and some expressions represent non-self-referentialvariable reassignment. Some final variable assignments to variables wereused in the final step of the calculation, and some final variableassignments to variables were not used in the final step of thecalculation. All calculations were executed from top to bottom, thenfrom left to right (fully computing each column before moving to thenext one).

All equations were verified for accuracy using a symbolic solver.

6.1.1 Common Subexpression Elimination of the Jacobian Matrix Elements

For each submatrix of the Jacobian, subexpressions were calculated, asshown in FIGS. 28 and 29 , as

$\begin{matrix}{a_{0} = {\overset{\_}{x}}_{ij}^{2}} \\{a_{1} = {\overset{\_}{y}}_{ij}^{2}} \\{a_{2} = {\overset{\_}{z}}_{ij}^{2}} \\{a_{7} = {a_{0} + a_{1} + a_{2}}} \\{a_{3} = {a_{7}\left( {{3a_{0}} - a_{7}} \right)}} \\{a_{4} = {a_{7}\left( {{3a_{1}} - a_{7}} \right)}} \\{a_{5} = {a_{7}\left( {{3a_{2}} - a_{7}} \right)}}\end{matrix}\begin{matrix}{a_{0} = {15a_{0}}} \\{a_{1} = {15a_{1}}} \\{a_{2} = {15a_{2}}} \\{b_{3} = {{\overset{\_}{x}}_{ij}{\overset{\_}{y}}_{ij}}} \\{a_{8} = {3a_{7}}} \\{b_{0} = {b_{3}a_{8}}} \\{b_{2} = {{\overset{\_}{z}}_{ij}a_{8}}} \\{b_{5} = {a_{0} - a_{8}}} \\{b_{7} = {a_{1} - a_{8}}} \\{b_{9} = {a_{2} - a_{8}}}\end{matrix}\begin{matrix}{a_{8} = {3a_{8}}} \\{a_{0} = {a_{0} - a_{8}}} \\{a_{1} = {a_{1} - a_{8}}} \\{a_{2} = {a_{2} - a_{8}}} \\{a_{6} = {{\overset{\_}{m}}_{j}a_{7}^{{- 7}/2}}}\end{matrix}$ $\begin{matrix}{b_{1} = {{\overset{\_}{y}}_{ij}b_{2}}} \\{b_{2} = {{\overset{\_}{x}}_{ij}b_{2}}} \\{b_{3} = {15b_{3}{\overset{\_}{z}}_{ij}}} \\{b_{4} = {{\overset{\_}{y}}_{ij}b_{5}}} \\{b_{5} = {{\overset{\_}{z}}_{ij}b_{5}}} \\{b_{6} = {{\overset{\_}{x}}_{ij}b_{7}}} \\{b_{7} = {{\overset{\_}{z}}_{ij}b_{7}}} \\{b_{8} = {{\overset{\_}{x}}_{ij}b_{9}}} \\{b_{9} = {{\overset{\_}{y}}_{ij}b_{9}}}\end{matrix}\begin{matrix}{c_{2} = \theta_{j}} \\{c_{3} = \phi_{j}} \\{c_{4} = {\sin c_{3}}} \\{c_{3} = {\cos c_{3}}} \\{c_{5} = {\sin c_{2}}} \\{c_{2} = {\cos c_{2}}} \\{c_{0} = {c_{8}c_{3}}} \\{c_{1} = {c_{5}c_{4}}} \\{c_{3} = {c_{3}c_{2}}} \\{c_{4} = {c_{4}c_{2}}}\end{matrix}$

giving variables equal to the expressions

$\begin{matrix}{a_{0} = {3\left( {{2{\overset{\_}{x}}_{ij}^{2}} - {3{\overset{\_}{y}}_{ij}^{2}} - {3{\overset{\_}{z}}_{ij}^{2}}} \right)}} \\{a_{1} = {3\left( {{{- 3}{\overset{\_}{x}}_{ij}^{2}} + {2{\overset{\_}{y}}_{ij}^{2}} - {3{\overset{\_}{z}}_{ij}^{2}}} \right)}} \\{a_{2} = {3\left( {{{- 3}{\overset{\_}{x}}_{ij}^{2}} - {3{\overset{\_}{y}}_{ij}^{2}} + {2{\overset{\_}{z}}_{ij}^{2}}} \right)}} \\{a_{3} = {\left( {{\overset{\_}{x}}_{i}^{2} + {\overset{\_}{y}}_{i}^{2} + {\overset{\_}{z}}_{i}^{2}} \right)\left( {{2{\overset{\_}{x}}_{ij}^{2}} - {\overset{\_}{y}}_{ij}^{2} - {\overset{\_}{z}}_{ij}^{2}} \right)}} \\{a_{4} = {\left( {{\overset{\_}{x}}_{i}^{2} + {\overset{\_}{y}}_{i}^{2} + {\overset{\_}{z}}_{i}^{2}} \right)\left( {{- {\overset{\_}{x}}_{ij}^{2}} + {2{\overset{\_}{y}}_{ij}^{2}} - {\overset{\_}{z}}_{ij}^{2}} \right)}} \\{a_{5} = {\left( {{\overset{\_}{x}}_{i}^{2} + {\overset{\_}{y}}_{i}^{2} + {\overset{\_}{z}}_{i}^{2}} \right)\left( {{- {\overset{\_}{x}}_{ij}^{2}} - {\overset{\_}{y}}_{ij}^{2} + {2{\overset{\_}{z}}_{ij}^{2}}} \right)}} \\{a_{6} = {{\overset{\_}{m}}_{j}\left( {{\overset{\_}{x}}_{ij}^{2} + {\overset{\_}{y}}_{ij}^{2} + {\overset{\_}{z}}_{ij}^{2}} \right)}^{{- 7}/2}}\end{matrix}\begin{matrix}{b_{0} = {3\left( {{\overset{\_}{x}}_{i}^{2} + {\overset{\_}{y}}_{i}^{2} + {\overset{\_}{z}}_{i}^{2}} \right){\overset{\_}{x}}_{ij}{\overset{\_}{y}}_{ij}}} \\{b_{1} = {3\left( {{\overset{\_}{x}}_{i}^{2} + {\overset{\_}{y}}_{i}^{2} + {\overset{\_}{z}}_{i}^{2}} \right){\overset{\_}{y}}_{ij}{\overset{\_}{z}}_{ij}}} \\{b_{2} = {3\left( {{\overset{\_}{x}}_{i}^{2} + {\overset{\_}{y}}_{i}^{2} + {\overset{\_}{z}}_{i}^{2}} \right){\overset{\_}{x}}_{ij}{\overset{\_}{z}}_{ij}}} \\{b_{3} = {15{\overset{\_}{x}}_{ij}{\overset{\_}{y}}_{ij}{\overset{\_}{z}}_{ij}}} \\{b_{4} = {3{{\overset{\_}{y}}_{ij}\left( {{4{\overset{\_}{x}}_{ij}^{2}} - {\overset{\_}{y}}_{ij}^{2} - {\overset{\_}{z}}_{ij}^{2}} \right)}}} \\{b_{5} = {3{{\overset{\_}{z}}_{ij}\left( {{4{\overset{\_}{x}}_{ij}^{2}} - {\overset{\_}{y}}_{ij}^{2} - {\overset{\_}{z}}_{ij}^{2}} \right)}}} \\{b_{6} = {3{{\overset{\_}{x}}_{ij}\left( {{- {\overset{\_}{x}}_{ij}^{2}} + {4{\overset{\_}{y}}_{ij}^{2}} - {\overset{\_}{z}}_{ij}^{2}} \right)}}} \\{b_{7} = {3{{\overset{\_}{z}}_{ij}\left( {{- {\overset{\_}{x}}_{ij}^{2}} + {4{\overset{\_}{y}}_{ij}^{2}} - {\overset{\_}{z}}_{ij}^{2}} \right)}}} \\{b_{8} = {3{{\overset{\_}{x}}_{ij}\left( {{- {\overset{\_}{x}}_{ij}^{2}} - {\overset{\_}{y}}_{ij}^{2} + {4{\overset{\_}{z}}_{ij}^{2}}} \right)}}} \\{b_{9} = {3{{\overset{\_}{y}}_{ij}\left( {{- {\overset{\_}{x}}_{ij}^{2}} - {\overset{\_}{y}}_{ij}^{2} + {4{\overset{\_}{z}}_{ij}^{2}}} \right)}}}\end{matrix}$ $\begin{matrix}{c_{0} = {\sin\theta_{j}\cos\phi_{j}}} \\{c_{1} = {\sin\theta_{j}\sin\phi_{j}}} \\{c_{2} = {\cos\theta_{j}}} \\{c_{3} = {\cos\theta_{j}\cos\phi_{j}}} \\{c_{4} = {\cos\theta_{j}\sin\phi_{j}}} \\{{c_{5} = {\sin\theta_{j}}},}\end{matrix}$

and yielding the submatrix

J_(ij) = a₀?. ?indicates text missing or illegible when filed

Assigning the remaining repeating elements of this submatrix to threeadditional local variables,

_(yz) =b ₃ c ₀ +b ₇ c ₁ +b _(g) c ₂

a _(zz) =b ₅ c ₀ +b ₃ c ₁ +b ₈ c ₂

a _(xy) =b ₄ c ₀ +b ₆ c ₁ +b ₃ c ₂,

gives Jacobian submatrices of the form

? ?indicates text missing or illegible when filed

noting the variable a₆ multiplied to all of the elements of the matrix.

When computing the full Jacobian matrix, it can be fastest to computethese submatrices column-wise first, noting that the orientationinformation contained in variables c₀ through c₅ is common through allsubmatrices corresponding to an individual magnet, and thus allowing thesines and cosines to be computed just once per magnet.

6.1.2 Common Subexpression Elimination of the Magnetic Field PredictionError Calculation

Now applying the same strategy as above to the magnetic field predictioncalculation, as shown in FIG. 30 , while as much as possible preservingsimilar notation as in the calculations above, subexpressions for themagnetic field prediction error were calculated as

$a_{0} = {\overset{\_}{x}}_{ij}^{2}$$a_{1} = {\overset{\_}{y}}_{ij}^{2}$$a_{2} = {\overset{\_}{z}}_{ij}^{2}$ a₇ = a₀ + a₁ + a₂$a_{8} = {{\overset{\_}{m}}_{j}{a_{7}}^{{- 5}/2}}$ c₂ = θ_(j) c₃ = ϕ_(j)c₄ = sin c₃ c₃ = cos c₃ c₅ = sin c₂ c₂ = cos c₂ c₀ = c₅c₃ c₁ = c₅c₄$d_{3} = {3\left( {{{\overset{\_}{x}}_{ij}c_{0}} + {{\overset{\_}{y}}_{ij}c_{1}} + {{\overset{\_}{z}}_{ij}c_{2}}} \right)}$$d_{0} = {{{\overset{\_}{x}}_{ij}d_{3}} - {a_{7}c_{0}}}$$d_{1} = {{{\overset{\_}{y}}_{ij}d_{3}} - {a_{7}c_{1}}}$${d_{2} = {{{\overset{\_}{z}}_{ij}d_{3}} - {a_{7}c_{2}}}},$

giving relevant local variables equal to the expressions

$a_{8} = \frac{{\overset{\_}{m}}_{j}}{\left( {{\overset{\_}{x}}_{ij}^{2} + {\overset{\_}{y}}_{ij}^{2} + {\overset{\_}{z}}_{ij}^{2}} \right)^{\frac{5}{2}}}$$d_{0} = {{3{{\overset{\_}{x}}_{ij}\left( {{{\overset{\_}{x}}_{ij}{\sin\left( \theta_{j} \right)}{\cos\left( \phi_{j} \right)}} + {{\overset{\_}{y}}_{ij}{\sin\left( \theta_{j} \right)}{\sin\left( \phi_{j} \right)}} + {{\overset{\_}{z}}_{ij}{\cos\left( \theta_{j} \right)}}} \right)}} - {\left( {{\overset{\_}{x}}_{ij}^{2} + {\overset{\_}{y}}_{ij}^{2} + {\overset{\_}{z}}_{ij}^{2}} \right){\sin\left( \theta_{j} \right)}{\cos\left( \phi_{j} \right)}}}$$d_{1} = {{3{{\overset{\_}{y}}_{ij}\left( {{{\overset{\_}{x}}_{ij}{\sin\left( \theta_{j} \right)}{\cos\left( \phi_{j} \right)}} + {{\overset{\_}{y}}_{ij}{\sin\left( \theta_{j} \right)}{\sin\left( \phi_{j} \right)}} + {{\overset{\_}{z}}_{ij}{\cos\left( \theta_{j} \right)}}} \right)}} - {\left( {{\overset{\_}{x}}_{ij}^{2} + {\overset{\_}{y}}_{ij}^{2} + {\overset{\_}{z}}_{ij}^{2}} \right){\sin\left( \theta_{j} \right)}{\sin\left( \phi_{j} \right)}}}$$d_{2} = {{3{{\overset{\_}{z}}_{ij}\left( {{{\overset{\_}{x}}_{ij}{\sin\left( \theta_{j} \right)}{\cos\left( \phi_{j} \right)}} + {{\overset{\_}{y}}_{ij}{\sin\left( \theta_{j} \right)}{\sin\left( \phi_{j} \right)}} + {{\overset{\_}{z}}_{ij}{\cos\left( \theta_{j} \right)}}} \right)}} - {\left( {{\overset{\_}{x}}_{ij}^{2} + {\overset{\_}{y}}_{ij}^{2} + {\overset{\_}{z}}_{ij}^{2}} \right){{\cos\left( \theta_{j} \right)}.}}}$

The magnetic field prediction error was then computed as

$\begin{matrix}{E_{ix} = {\left( {G_{x} + {\sum\limits_{j^{\prime} = 0}^{j^{\prime} = M}{d_{0}a_{8}}}} \right) - {\overset{\sim}{B}}_{ix}}} & \left( {6.2a} \right)\end{matrix}$ $\begin{matrix}{E_{iy} = {\left( {G_{y} + {\sum\limits_{j^{\prime} = 0}^{j^{\prime} = M}{d_{1}a_{8}}}} \right) - {\overset{\sim}{B}}_{iy}}} & \left( {6.2b} \right)\end{matrix}$ $\begin{matrix}{E_{iz} = {\left( {G_{z} + {\sum\limits^{j^{\prime} = M}{d_{2}a_{8}}}} \right) - {\overset{\sim}{B}}_{iz}}} & \left( {6.2c} \right)\end{matrix}$

Following a similar strategy to that above, it can be best to computethe prediction of the magnetic field contributions from each magnetfirst, to avoid recomputing sines and cosines. Magnetic fieldcontributions from each magnet were summed to vector components as theywere computed, and then, afterwards, the disturbance field and measuredfield were added and subtracted, respectively.

6.2 Practical Notes

Common notation was deliberately used between the subexpressions for themagnetic field prediction vector and the subexpressions for the Jacobianmatrix. As indicated by this common notation, further commonsubexpression elimination can be performed by combining the magneticfield prediction vector and Jacobian matrix calculations, furtherreducing the latency of the optimization. However, for simplicity, inthis work, it was chosen not to combine these two calculations, becausethe optimization algorithm employed, the Levenberg-Marquardt algorithm,does not request the optimization function and its derivative the samenumber of times. To combine the calculations, the source code of thisalgorithm can be altered to ensure that the magnetic field predictionvalues are not overwritten whenever only the derivative is requested.Use with other optimization algorithms or other hardware may makecombining these two calculations more clearly advantageous.

In addition, parallelization was found to be inefficient with thehardware used in this work, but the computation of the magnetic fieldprediction error Jacobian matrix is a parallel process. With specializedhardware, the Jacobian matrix can be calculated using up to MN differentprocessing streams, and the part of the magnetic field predictioncalculation before the summation step can also be parallelized alongwith the Jacobian matrix. Further, an inspection of the commonsubexpression elimination flow diagrams (see FIGS. 28-30 ) reveals thatthese calculations are well-suited for implementation on an ASIC(application-specific integrated circuit) or an FPGA (field-programmablegate array) using hardware description language, on which fullparallelization without delays from overhead can be possible.

A significant part of the time delay in these calculations is due to thesine and cosine operations, so any sizeable advance in the speed ofthese operations can have a non-negligible effect on the speed of thetracking algorithm. Noting that the sine and cosine operations in thesesubexpressions are both performed on both θ_(j) and φ_(j), and furthernoting that there is some redundancy of information between sine andcosine, the computational latency was able to be

further reduced by using the sincos( ) function [21]. This modificationwas implemented after the results of the Taylor Reference, but beforethe experiments described in Example 1 herein. There are other methodsfor combining the sine and cosine operations, but the use is againhardware-specific. For example, an implementation of CORDIC tested withour algorithm was found to be slow on macOS, but CORDIC may be the bestchoice on another platform.

Depending upon the particular optimization implementation, the Jacobianmatrix can be transposed or vectorized from its current structure beforeusage with an optimization algorithm. With this current implementation,it was found to be fastest to store the Jacobian directly as a vectorinstead of using multiple dimensions, allowing storage and access withfewer operations. The final step of the common subexpression eliminationtechniques discussed above can be implemented in a manner which directlyproduces this vectorization of the Jacobian using a sparse matrixmultiplied by a vector containing the orientation information, butwithout extremely specialized hardware, this would likely not be asefficient as computing each of the elements as discussed above,unnecessarily requiring additional memory and increasing the number ofcommon subexpressions required.

While every effort was made to minimize the number of local variables ineach function, in the case of the products boco and boci, where theseoperations are only repeated when tracking the dipole strength, theseduplicate subexpressions were allowed to remain. This decision was justfor simplicity, noting that low latencies were only required whentracking magnetic dipoles of known size. Similarly, when trackingmagnets with reduced dimensionality (e.g., position tracking on asurface with fixed orientation), some of the intermediate variables usedabove need not be computed, and the particular common subexpressionelimination method can be adapted to the particular degrees of freedomused.

6.3 Improvement Via Newton's Method

It is worthy of considerable focus to investigate the use ofsecond-order optimization using Newton's Method of optimization. InChapter 2 and Chapter 3, we used the Levenberg-Marquardt algorithm tominimize the sum of the squared magnetic field prediction error,

$\begin{matrix}{{S = {\sum\limits_{i = 1}^{N}{\sum\limits_{\underset{\{{x,y,z}\}}{*=}}^{}E_{i*}^{2}}}},} & (6.3)\end{matrix}$

subject to the tracking parameters

p=(x ₁ ,y ₁ ,z ₁,θ₁,ϕ₁ ,{umlaut over (m)} ₁ , . . . ,x _(M) ,y _(M) ,z_(M),θ_(M),ϕ_(M) ,{umlaut over (m)} _(M) ,G _(x) ,G _(y) ,G _(z)).

The first order partial derivative of S with respect to the jth trackingparameter is given by

$\begin{matrix}{{\left( {\nabla S} \right)_{j} = {2{\sum\limits_{i = 1}^{N}{\sum\limits_{\underset{\{{x,y,z}\}}{*=}}^{}{E_{i*}\frac{\partial E_{i*}}{\partial p_{j}}}}}}},} & (6.4)\end{matrix}$

which is the (doubled) dot product of the jth column of the Jacobianmatrix of E with E, where E=(E_(1x),E_(1y),E_(1z), . . .,E_(Nx),E_(Ny),E_(Nz))|. More succinctly, these partial derivatives canbe expressed as the vector

∇S=2J^(T)E.  (6.5)

This vector ∇S was the gradient of the optimization used for magnettracking with the LevenbergMarquardt algorithm.

To use Newton's Method, the second order partial derivatives of S,corresponding to the matrix elements of the Hessian of S, is alsoprovided

$\begin{matrix}{H_{jk} = {2{\sum\limits_{i = 1}^{N}{\sum\limits_{\underset{\{{x,y,z}\}}{*=}}^{}{\left( {{\frac{\partial E_{i*}}{\partial p_{j}}\frac{\partial E_{i*}}{\partial p_{k}}} + {E_{i*}\frac{\partial^{2}E_{i*}}{{\partial p_{j}}{\partial p_{k}}}}} \right).}}}}} & (6.6)\end{matrix}$

The summed first term of this equation is simply the dot product of thejth column of the Jacobian matrix of E with the kth column of theJacobian matrix of E. The summed second term is an element of the tensorcontraction of the Hessian of E with E. In other words, if the elementsof the Hessian of E are given by

$\begin{matrix}{{{\overset{\prime}{H}}_{i*{jk}} = \frac{\partial^{2}E_{i*}}{{\partial p_{j}}{\partial p_{k}}}},} & (6.7)\end{matrix}$

then (6.6) can be rewritten as

$\begin{matrix}{{H_{jk} = {2{\sum\limits_{i = 1}^{N}{\sum\limits_{\underset{\{{x,y,z}\}}{*=}}^{}\left( {{\frac{\partial E_{i*}}{\partial p_{j}}\frac{\partial E_{i*}}{\partial p_{k}}} + {E_{i*}{\overset{\prime}{H}}_{i*{jk}}}} \right)}}}},} & (6.8)\end{matrix}$

and the Hessian of S is given by the matrix equation

H=2(J ^(T) J+EH′),  (6.9)

where H′ is the Hessian of E, and EH′ represents the tensor contractionof the Hessian of E with E. Representing this tensor contraction asH=EH′ and defining the “half-Hessian” of S as H Δ ⅓ H gives

H ⁻ =J ^(T) J+H.  (6.10)

Borrowing scale invariance from the Levenberg-Marquardt algorithm, theequation for the vector step size δ can then be represented as

(H ⁻+λdiag(H ⁻))δ=J ^(T) E.  (6.11)

Upon solving (6.11) for δ, the tracking parameters can then beiteratively updated as

p _(t+1) =p _(t)−δ  (6.12)

Newton's method, when combined with the scale invariance parameter fromthe Levenberg-Marquardt algorithm as in (6.11), can exhibit increasedstability in the magnet tracking problem. Further, it is anticipatedthat the additional computational delay required to calculate theHessian elements can be compensated for by a reduced step count neededto reach convergence.

Common subexpression elimination (and parallelization, with applicablehardware) can substantially accelerate the calculation of the Hessianmatrix, making the solving of the matrix equation (6.11) the dominantdriver of time-delay. The use of common subexpression elimination can beapplied simultaneously to the calculation of the Hessian matrix,Jacobian matrix, and magnetic field predictions. Further, the simplicityof executing the matrix algebra for Newton's method can open thepotential to implement the entire algorithm on an FPGA. As such, thetime delay of the tracking algorithm can be substantially furtherreduced.

7. Further Advances

To address the high magnet tracking time delays in state-of-the-artmultiple-magnet tracking, a high-speed multiple-magnetic-target trackingalgorithm was developed, with latencies low enough to enable reflexiveprosthesis control. Then, finding limitations to the use of magnettracking in portable applications, a disturbance rejection algorithm wasdeveloped which can be implemented entirely in software. To address theneed for a framework to calibrate multiple magnetic field sensors atonce, previous ellipsoid fitting techniques were extended and improvedto adapt them specifically to multi-sensor calibration. Lastly, seeing aneed for a minimally-invasive technology to accurately track musclelengths and speeds in real time, a new technology, magnetomicrometry wasdeveloped to meet this need, and the technology was validated bytracking muscle lengths in an animal model.

A framework for the next generation of magnetomicrometry is provided,which comes with a host of additional improvements. The orientation andposition calibration discussed in Section 4 herein and the improvedmagnetic dipole strength measurement discussed in Sections 3 and 4 canenable new sensor geometries and improved accuracy. For example,magnetic field sensors can be placed on a curved array (see FIG. 34 ),stacked in multiple layers, or both, without a need to precisionassemble or optically measure the positions and orientations of thecomponents. Assembly of a curved array, such as shown in FIG. 34 , can,for instance, be performed by fabricating the circuit on a flexibleboard and using an adhesive to affix it to a rigid non-planar geometry.For example, epoxy can be used to attach flexible magnetic field sensingcircuitry flush with an outer surface of a prosthetic socket, ensuringthat the magnetic field sensors are placed as proximal to the magneticbeads as possible. This flexible magnetic field sensing circuitry cantake a form similar to the multiple sensor grids fabricated for theexperiments described in Section 1 and Example 1 herein, or it can takea form similar to an LED tape light strip, allowing it to becontinuously wrapped around the prosthetic socket.

An example of a non-planar sensing array 1100 is shown in FIG. 34 . Asillustrated, the non-planar sensing array is curved, but othergeometries are possible, including irregularly curved and angleddevices. The array 1100 includes a plurality of sensing elements 1101and can be positioned against anatomy, or in or on a prosthetic devicethat encloses an anatomy (not shown for clarity), in which magnetictargets 101 are implanted.

A method of assembling a non-planar sensing array includes fabricating aplurality of sensors (e.g., magnetic field sensors) on a flexiblecircuit board and affixing the flexible circuit board to a non-planarand rigid substrate (e.g., a prosthetic socket).

7.1 Sensing of Systems of Permanent Magnets Wherein the RelativePositions of Multiple Sensing Arrays are Tracked via InertialMeasurement Units

One advantage of using multiple magnetic beads within a single tissue(e.g., muscle) is that this sensing modality is not affected bymovements of the sensing array relative to the muscle. When multiplesensing arrays are used, these sensing arrays can be rigidly fixed inposition and orientation relative to one another to ensure that movementof the arrays does not result in sensing noise. For example, themultiple arrays can be fixed to a carbon fiber prosthetic socket. Theremay arise situations, however, in which it is best to not have a rigidelement holding the arrays together. For example, for an amputee with anosseo-integrated prosthesis, a socket is not needed. In this instance,it may be better to have the sensing arrays attached independently.Methods and systems providing for independent, flexible attachment ofmultiple sensing arrays to track systems of magnetic beads are provided.In such systems and methods, not all magnetic beads need be sensed byall sensing arrays in order for tracking to occur of each bead relativeto the system of beads.

Such methods and systems provide for the use of multiple sensor arraysand in which the arrays are tracked relative to one another.

An example system in shown in FIG. 31 . In this example, more than onesensing array is included (e.g., sensing arrays 1208 a, 1208 b, eachwith a plurality of sensing elements 1201 configured to detect a stateof at least one magnetic target at a tissue). Each sensing array 1208 a,1208 b includes a position sensor 1210 (e.g., an inertial measurementunit, accelerometer, etc.) to monitor a position and orientation of eachsensing array (hereafter referred to as the pose of the sensing array)relative to the at least one other array. A controller 1230 in a wiredor wireless arrangement with the sensing arrays can be configured toreceive position data of the sensing arrays.

The pose of each array can be compared against the pose of the otherarrays. For example, one of the arrays can be designated as thereference array, with a global coordinate system defined to be that seenby the reference array. The pose of all other sensing arrays can then becalculated with respect to the reference array. Regardless of how thecoordinate system is defined, all sensing arrays of the system can bebrought into a common coordinate system. More particularly, thepositions and orientations of each of the magnetic field sensors on eachsensing array are then calculated in this new common coordinate systemusing the calculated relative pose of the sensing arrays, and themagnetic beads are then tracked with respect to the common coordinatesystem.

An example method is shown in FIG. 32 . The method 1300 includes thesensor arrays beginning at a known position and orientation relative toone another (1302). Translation acceleration and angular velocity assensed by an inertial measurement unit (IMU) of each array can be usedto monitor the relative poses of the sensor arrays (1304) and, inparallel, each array can track at least a subset of the magnetic targetsincluded in or at an anatomy (1306). The relative poses of the sensorarrays and tracking parameters for at least a subset of magnets trackedby each array are combined to determine relative positions of allmagnets in the system relative to one another (1308).

The sensor arrays can be connected to one another in a manner thatrestricts their relative degrees of freedom and enables the remainingrelative degrees of freedom to be sensed with greater simplicity. Forexample, one or more hinges can be used for two or more sensor arrays,with one or more potentiometers to sense the angle between the sensorarrays. In another example, one or more linear bearings can be used toreduce the relative movement of the arrays to a single degree offreedom, and the relative positions can be measured via a distancesensor, such as an infrared sensor. In either of the aboveimplementations, a flex sensor, IMU, or other sensing strategy can beused to determine the relative poses of the sensor arrays.

7.4 Preferred Sensor Magnetic Bead Tracking

When tracking magnetic beads, there can be a tradeoff between accuracyand tracking speed. From a hardware perspective, this tradeoff can bedriven by the number of sensors that are used when tracking and thenumber of magnets that are tracked simultaneously by a shared set ofsensors. This section describes an example system and method in whichthis tradeoff can be overcome by tracking multiple magnets for differentsets of sensors that are themselves monitored with respect to theirrelative locations. The following describes a method whereby thistradeoff can be overcome by adjusting how the magnetic field sensors areused in the tracking.

A method is described for the strategic use of subsets of sensor datawhich can result in improved tracking accuracy. By using highly relevantsubsets of the sensor data for tracking, the tracking accuracy and thesize of the tracking volume can be improved upon with reduced trackingrequirements.

For example, with the system 1200, the controller 1230 can further beconfigured to determine a distance and difference in orientation betweeneach of the at least two sensing arrays and at least one magnetictarget. Optionally, at least one magnetic target can be associated withone of the at least two sensor arrays based upon the determineddistances and orientations.

7.4.1 Nearest Sensor Ranking for Magnetic Bead Tracking, Via TrackedMagnetic Bead Location

In this example, the tracked location of the magnetic beads is used torank the relevance of the sensors. For each bead at each timepoint, thelocation of the bead as tracked at the previous timepoint is used (orthe previous position, velocity, orientation, angular velocity, and/oraccelerations of the bead can be used to predict the current timepointparameters), the distance between the bead and the sensor is calculated,and all sensors are ranked according to their distance from the trackedlocation of the bead, with shorter distances imparting higher ranking.For efficiency, the square of the distance may instead be used toperform this ranking. To rank the sensors when multiple beads are used,the ranking can be performed for each of the beads independently and thetwo rankings can be merged. There are many different ways that theserankings can be merged, but one example of this merging iterates betweenthe two lists in compiling a new list while ensuring that eachindependent list makes an equal (or off-by-one) number of contributionsto the merged list. To track the magnets using a subset of the set ofsensors, the subset is chosen giving priority to higher ranked sensorsfirst.

7.4.2 Strongest Sensor Ranking for Magnetic Bead Tracking, Via EstimatedMagnetic Bead Magnetic Field Contribution at Each Sensor

While the example described in Section 7.4.1 can provide a simple,efficient rule of thumb for selecting a subset of sensors, it may notyield the highest accuracy. Instead, the sensors that are the mostrelevant with respect to each magnetic bead can be considered to be thesensors which measure the greatest magnetic fields. The example methodand system described herein addresses this relevance by ranking thesensors based upon the magnetic field contribution of each of themagnetic beads. In this example, for each bead at each timepoint, thelocation of the bead as tracked at the previous timepoint is used, andthe predicted magnetic field is calculated along each axis of eachsensor. For each magnetic bead, the sensors are ranked according to thestrength (magnitude) of the magnetic field as measured at the sensors.This ranking can be performed using the magnitude of the vector magneticfield predicted from all three axes at once, or it can be performedranking all axes of all sensors independently (resulting in trackingwhich does not necessarily include all three components of each sensorused in the tracking). When multiple magnetic beads are tracked, theindependent rankings can be combined as described in Section 7.4.1.Again, in this instance, the subset of the set of sensors chosen can beselected based upon priority, with higher ranked sensors being selectedfirst.

A flowchart is shown in FIG. 33 illustrating an example method 1400 ofranking for magnetic bead tracking. The method includes, for eachtracked magnetic bead (1402), predicting the current position andorientation based on the previous position, orientation, andtranslational and angular velocity of the bead (1404). The methodfurther includes calculating a predicted magnetic field contributionfrom the considered bead at the predicted position (1406) and rankingthe sensors of the system based on the magnitudes of the vectorpredicted magnetic fields at each sensor (1408). For a number ofiterations smaller than the number of sensors, and starting with thefirst bead (1410), the highest-ranking sensor corresponding to the beadis added to a global ranking list and the sensor is removed from thelocal ranking list (1412). With the beads cycled through in a consistentorder, a next bead is used for a next iteration (1414).

7.5 Further Considerations

In addition to improved geometry, improved tracking methods such as theuse of the Hessian to implement Newton's method (as discussed in Section6.3 herein) can allow for further increased tracking stability, and canfurther reduce magnet tracking time delay. The use of specializedhardware, such as a hybrid or full FPGA approach, can substantiallyreduce the time delay required for magnet tracking. Increased trackingstability and reduced tracking time delay further allow tracking moremagnets simultaneously in real time, allowing additional flexibility insensor positioning when tracking many magnets at once.

To further maintain both accuracy and time delay while extendingtracking range over a larger volume, an algorithm that makes use of onlythe nearest sensors to each magnet can be used. This could, forinstance, first use all sensors to determine the initial magnetlocations, and then use the tracked magnet positions to sort the sensorsaccording to magnet proximity. Alternatively, the signal-to-noise ratiofor each sensor can be ranked based on each magnet's predicted magneticfield.

In theory, an instability should exist when the magnet is aligned alongits reference axis. Though no practical issues were encountered aroundvertical orientation in this work, if instability along the referenceaxis is determined to be non-negligible, the use of quaternions may notbe sufficient to solve this issue due to the symmetry of the magnet.However, using an alternative reference axis whenever the magnet isoriented near the primary reference axis can solve this issue, if itexists. The use of the Hessian matrix via Newton's method ofoptimization can be sufficient to remove any potential of theinstability.

No clear method for estimating the error of magnet tracking on-the-flyhas been proposed, but it is possible to use the final Jacobian matrixof each tracking step to estimate the error of the current magnetposition and orientation estimate using propagation of uncertainty.Because the Jacobian matrix gives the change in magnetic field errorwith respect to position and orientation parameters, each element of theJacobian matrix can be inverted to give the error in position andorientation parameters that would be induced by an error in the magneticfield measurement. Assuming that the magnetic field error has a known,normal distribution, these inverted Jacobian matrix elements can then beweighted by the standard deviation of this normal distribution andaveraged to compute an estimate of the error. If more precision in theerror estimate is needed (including in the asymmetry of the errorestimate about the estimated pose), the Hessian elements can be used tofurther refine this error estimate using the inverse of the secondderivative elements. These error estimates can then be used as part of aKalman filter to, in turn, increase the tracking accuracy.

All sensors are inherently noisy. When magnets are tracked using anoptimization algorithm, the algorithm performs the equivalent functionof smoothing out this noise, combining the data from multiple sensors todetermine a single estimate for each parameter. Assuming that themagnetic field sensor noise is normally distributed (that the noise ischiefly Johnson and flicker noise), and that the error in a singleposition or orientation parameter that would be caused by one individualcomponent of a three-axis sensor in isolation is linearly related to theerror in sensing, it would then be expected that there exists acorresponding error distribution for each tracking parameter, comprisingerrors which correspond to each component of each three axis sensor.Though the optimization used for magnet tracking approaches the magnetparameter estimation holistically rather than individually, thesetheoretical individual tracking parameter error distributions can beteased out and one can then calculate their standard deviations and usethese standard deviations as a proxy for error estimation of each of theparameters during each tracking step.

To determine how these individual error estimates can be calculated, wenote that the Jacobian matrix in our magnet tracking algorithm providesthe derivative of the magnetic field prediction, and that at the end ofeach tracking step we have available the estimated magnet parameters(e.g., x₁) and both the predicted magnetic field (e.g., B_(1(1)y),referring to the predicted magnetic field contribution from the firstmagnet on the y-axis component of the first sensor) and measuredmagnetic field (e.g. {tilde over (B)}_(1y)) values for each component ofeach sensor. For a single component of a sensor and a single magnetparameter, then, when tracking a single magnet, we have for example, inpoint-slope form the equation

$\begin{matrix}{{B_{1{(1)}y} - {\overset{\sim}{B}}_{1y}} = {\frac{\partial B_{1y}}{\partial x_{1}}\left( {x_{1} - {\overset{\sim}{x}}_{1{({1y})}}} \right)}} & (7.1)\end{matrix}$

where x^({tilde over ( )}) _(1(1y)) is the magnet parameter as correctedusing the measured magnetic field. If we then let

e _(x1(1y)) =x ₁ −x ^({tilde over ( )}) _(1(1y))  (7.2)

represent the estimated error in estimate of xi corresponding to themeasured magnetic field B^({tilde over ( )}) _(1y), we have

$\begin{matrix}{e_{x1{({1y})}} = \frac{B_{1{(1)}y} - {\overset{\sim}{B}}_{1y}}{\frac{\partial B_{1y}}{\partial x_{1}}}} & (7.3)\end{matrix}$

and in generalizing these particular parameters to the context ofmultiple magnet tracking, we can subtract the predicted magnetic fieldcontribution from all magnets other than the first magnet to retainvalidity of the equation as

$\begin{matrix}{e_{x1{({1y})}} = \frac{B_{1{(1)}y} - \left( {{\overset{\sim}{B}}_{1y} - \left( {B_{1y} - B_{1{(1)}y}} \right)} \right)}{\frac{\partial B_{1y}}{\partial x_{1}}}} & (7.4)\end{matrix}$ $\begin{matrix}{= \frac{B_{1y} - {\overset{\sim}{B}}_{1y}}{\frac{\partial B_{1y}}{\partial x_{1}}}} & (7.5)\end{matrix}$ $\begin{matrix}{= \frac{E_{1y}}{\left( \frac{\partial B_{1y}}{\partial x_{1}} \right)}} & (7.6)\end{matrix}$

This gives us a distribution of errors of our singled-out parametercorresponding to all of the sensor components, which distribution we canthen use to better understand the accuracy of our parameter estimation.For instance, we can then calculation the standard deviation of thedistribution of the estimated errors in the x₁ estimate as

$\begin{matrix}{e_{x1_{UNWEIGHTED}} = \sqrt{\frac{e_{x1{({1x})}}^{2} + e_{x1{({1y})}}^{2} + e_{x1{({1z})}}^{2} + \cdots + e_{x1{({Nx})}}^{2} + e_{x1{({Ny})}}^{2} + e_{x1{({Nz})}}^{2}}{3N}}} & (7.7)\end{matrix}$

A vector of all of the standard deviations in this unweighted form canbe calculated via the generalized equation

$\begin{matrix}{e_{UNWEIGHTED} = {\left( {\frac{1}{3N}\left( {J{^\circ}^{- 2}} \right)^{T}E{^\circ}^{2}} \right){^\circ}^{\frac{1}{2}}}} & (7.8)\end{matrix}$

where

represents element-wise (Hadamard) operations on the matrices. Toaccount for the skewing that occurs in the values that are far away, theweighting

$\begin{matrix}{w_{ij} = \frac{r_{ij}^{- 4}}{3{\sum}_{i = 1}^{N}r_{ij}^{- 4}}} & (7.9)\end{matrix}$

can be applied to each submatrix corresponding to i and j of theelement-wise-squared-inverse Jacobian matrix instead of division by 3N.

In summary, to determine an estimate of the error for each of thetracking parameters, at the end of each tracking step, the Jacobian canbe inverted element-wise and squared element-wise, weighted by theinverse fourth power of the distance from the magnets to the sensors toaccount for skewing, then multiplied by the element-wise squaredmagnetic field prediction error column-vector. The element-wise squareroot of the resulting column-vector should then be used in estimatingthe tracking parameter errors.

Even without these improvements, magnetomicrometry is well-positioned tohave a tremendous impact in many fields. However, these next-generationmodifications will serve to increase the robustness, speed, and accuracyof this tool, enabling it to have an even farther-reaching effect.

Additional descriptions of systems and methods relating to tracking ofmagnetic objects can be found in WO2019/074950, “Method forNeuromechanical and Neuroelectromagnetic Mitigation of Limb Pathology;”the teachings of which are incorporated herein in their entirety.

Additional descriptions of systems and methods relating to tracking ofmultiple targets, such as permanent magnets, can be found in thefollowing publication: Cameron R Taylor, Haley G Abramson, and Hugh MHerr. Low-latency tracking of multiple permanent magnets. IEEE SensorsJournal, 19(23):11458-11468, 2019; the teachings of which areincorporated herein in their entirety.

EXEMPLIFICATION Example 1 Minimally-Invasive Muscle Tracking UsingPermanent Magnets

An in-vivo turkey model with a prototype magnetomicrometry system wasinvestigated to address biocompatibility, verify in-vivo trackingaccuracy, and ensure long-term resistance of the implants to migration.

Methods

All animal experiments were approved by the Institutional Animal Careand Use Committees at Brown University and the Massachusetts Instituteof Technology.

For surgical implantation of the magnetic beads, four turkeys wereplaced on anesthesia under 3-4% isofluorane, then plucked on both legsand scrubbed with Povidone-Iodine 10% (PVP-I). Turkeys were thentransferred to a hot pad, intubated and actively ventilated, whilemonitoring pulse oximetry, breathing rate, and body temperature. Thealternate leg was then wrapped in a sterile covering, the perimeter ofthe surgical field was draped with sterile cloth, and the leg wassterilized with PVP-I. At each insertion site (the distal and proximalends of the gastrocnemius and iliotibialis cranialis muscles), a 16gauge needle and a thin pair of unopened surgical scissors wereconsecutively used to make an insertion channel smaller than thediameter of the magnet. The magnet was then press-fit into the end of asterile hollow plastic tube, dipped in saline, and inserted into thechannel using depth markings on the plastic tube for reference. Asterile wooden rod (longer than the plastic tube) was then guided fullyinto the bore of the plastic tube and used to gently, but firmly, holdthe magnet in place while removing the plastic tube from the muscle. Thewooden rod was then removed, and nonmagnetic forceps were used to suturethe muscle closed at the insertion site using 6-0 nonabsorbable silksuture for the muscle and 4-0 Vicryl absorbable suture for the skin,after which Tegaderm (3M) transparent film dressing was applied to theskin around the insertion site.

For biocompatibility, all magnets (3-millimeter-diameter N35neodymium-iron-boron spherical magnets, initially coated in nickel) werecoated in parylene C (BJA Magnetics). Each magnet's strength was thenmeasured and recorded, and the magnet was rinsed in 70% ethanol byvolume (in distilled water) followed by three rinses with distilledde-ionized water. Each magnet was then placed in a labeled, ventedcontainer and sterilized using ethylene oxide, after which they wereallowed 48 hours to degas before surgical implantation.

During surgical implantation, magnet pairs were inserted, with the aidof a sterile ruler, at various separation distances betweenapproximately 20 and 70 millimeters. Immediately after surgicalimplantation and at time intervals (multiple weeks) after theimplantation, computed tomography (CT) scans (Animage Fidex VeterinaryCT Scanner) were used to monitor the distances between the beads.Turkeys were placed on anesthesia under 3-4% isofluorane, and for eachleg, the turkey lay prone with the leg of interest flush with, centeredon, and parallel to the scanning table, with the foot positioned ascranial and medial to the body as possible. Each leg was scannedseparately to simplify positioning in the scanner and reduce thepossibility of needing to repeat scans. In each CT scan, a referenceobject (an acrylic bar with magnets press-fit into two measured,pre-drilled holes) was included to ensure consistency in scale. Amedical image viewer (Horos) was used to determine the three-dimensionalpositions of the magnetic beads in each muscle, and these positions wereused to calculate the magnetic bead separation distances.

To track the magnets, two custom sensing boards, each with 48 LIS3MDLmagnetic field sensors (STMicroelectronics) in a 6×8 grid spaced at 4.83millimeters, were held together by a 3d-printed fixture (Connex 500,Stratasys) at 60 millimeters apart from board center to board center,forming a single, 96-magneticfield-sensor array (see FIGS. 3A-3B). Nylonnuts and bolts (McMaster-Carr) were used to secure the sensing boards tothe fixture. A custom adapter board was used to connect a Teensy 3.6microcontroller (PJRC) to the sensing boards using flexible flat cables(Molex), and on-board 4-to-16 line decoders (74HC154BQ, Nexperia) wereused to individually enable magnetic field sensors for SPI communication(10 MHz clock).

The magnetic field at each of the sensors was measured with a samplingrate of 300 Hz. To minimize onboard magnetic field distortion, allcapacitors used (Vishay) were MRI-safe. The tracking algorithm was runin real-time on a Macbook Air (13-inch, Early 2014) with 8 GB of RAM andan Intel i7 CPU running at 1.7 GHz.

To validate magnetomicrometry (see FIG. 4 ), for each of the legs ofeach of the four turkeys, the 96-magnetic-field-sensor array 250 a wassecured with veterinary tape to the outside of the turkey's leg over themagnetic bead pair 270 in the gastrocnemius muscle. With the turkeyanesthetized (3-4% isofluorane), an electric motor 272 (AuroraScientific 310B-LR) was used to apply a mechanical frequency sweep tothe turkey's ankle (10-second exponential chirp from 0.7 to 7 Hz), witha spring 274 (surgical tubing) providing an opposing force. Throughoutthis frequency sweep of the ankle (and thus, of the passively-cycledgastrocnemius muscle), the magnetic field sensor array was used, asdescribed in the Taylor Reference, to track the length of thegastrocnemius muscle using the distance between the magnetic beads inreal time. For comparison, the distance between the magnetic beads,which are radio-opaque, was also simultaneously monitored viafluoromicrometry (stereo X-ray videofluoroscopy), with the X-ray sources282 above the turkey and the image intensifiers 280 positioned below.All fluoromicrometry data was post-processed in XMALab, wheneverpossible automating the processing using 25% “threshold offset inpercent,” manually performing tracking when reprojection error exceededone pixel, and without performing any temporal filtering to smooth thedata. Time syncing was used to perform initial alignment of themagnetomicrometry and fluoromicrometry curves, but due to inconsistencyin the time sync signal from the X-ray system, optimization was used tofine-tune the temporal alignment of the data.

To investigate static tracking between magnetomicrometry andfluoromicrometry, for evaluating compatibility between the twotechnologies as well as sensing depth for the magnetomicrometry, twomagnets were placed into a 1×10 LEGO plate at various known distancesapart from one another, and the magnetic field sensor array was placedabove the magnets at various sensing heights.

To then further investigate the presence of translational offsetsbetween the sensor arrays, and to further investigate sensing depth, thetwo magnets were again placed into the 1×10 LEGO plate at 24 millimetersapart from one another, with a single magnetic field sensor board abovethe magnets at a finer interval of sensing heights.

Turkeys were labeled with identification cuffs by number and color,Green #5 (Turkey A), BlueGreen #2 (Turkey B), Yellow #1 (Turkey C), andBlue #6 (Turkey D).

Results: Magnetomicrometry

The results of the in-vivo turkey magnetomicrometry studies are shown inFIGS. 5-9 , with FIGS. 5-8 showing the distance between the magnets overtime throughout the frequency sweep for each trial of each leg of eachturkey. The experiments corresponding to Turkey A (FIGS. 5A-D) representfrequency sweeps on the right leg only, because the magnetic bead pairin the left gastrocnemius of Turkey A migrated together, as describedbelow with respect to long-term stability of the implants againstmigration. In the remaining turkeys, B through D (FIGS. 6A-F, 7A-F, and8A-F), plots show the three frequency sweeps corresponding to each leg.In all plots, magnetomicrometry is shown against fluoromicrometry, aswell as the absolute difference between them. In FIG. 9 , thedifferences between magnetomicrometry and fluoromicrometry are given, inhistogram form, for all of the results shown in FIGS. 5-8 . Inspectionof these figures shows a consistent offset between magnetomicrometry andfluoromicrometry for each set of trials for a given leg. This offset ismost clearly demonstrated in the trials for Turkeys C and D. It isexpected that this issue is due to misalignment of the magnetic fieldsensing boards, because they were attached to one another using a3D-printed fixture and plastic screws. It is thus expected that theseconsistent offsets seen can be corrected with calibration orhigh-precision fixturing. As shown by the table in FIG. 9 , the standarddeviation of the differences is consistently below 100 micrometers,suggesting that with correction for positioning, the system is capableof at least this precision.

To further compare magnetomicrometry and fluoromicrometry in acontrolled, static setting, ensuring that the two technologies did notinterfere with one another and investigating the need for proximity ofthe sensors to the magnets, the magnetic field sensor array was set atvarious heights and the magnetic beads were separated to variousdistances, and both technologies were used to simultaneously measure thedistances between the magnets. FIG. 10 shows these results, where themagnetic beads were placed at distances of approximately 24, 40, 56, and72 mm from one another and the sensor was placed at various heightsabove the magnetic beads. Note that, with the magnetic field sensorsclose to the magnetic beads, magnetomicrometry has a measurementstandard deviation that is significantly lower than fluoromicrometry,but as the magnetic field sensors are distanced from the beads, themagnetomicrometry measurement standard deviation approaches and thenbecomes larger than that of the fluoromicrometry measurement standarddeviation. Also note that, again, with the exception of the 40 mmmagnetic bead separation measured with a 30.5 mm sensor proximity, thereis a consistent offset between the magnetomicrometry andfluoromicrometry measurements, again suggesting the need for precisionfixturing or sensor position calibration.

Visualizing this same data temporally (see FIG. 11 ), we find that thefluoromicrometry does not interfere with the magnetomicrometrymeasurements (the magnetomicrometry measurements are consistent before,during, and after the X-ray irradiation).

Finally, to further investigate the offset that was encountered betweenthe magnetomicrometry and fluoromicrometry measurements, and to betterunderstand accuracy at various magnetic field sensor proximities, thestatic experiment was repeated using just one single board (withoutfluoromicrometry) above magnetic beads spaced approximately 24 mm apartin the 1×10 LEGO plate (see FIG. 12 for the results). Using a singleboard ensures that the positions of the magnetic field sensors are allwell-known relative to one another. In this experiment, separationdistance measurements were found to have high-accuracy at near proximityand sub-millimeter accuracy continuing to a sensing distance of justover 30 mm, followed by a pronounced increase in error thereafter.

Results: Long-Term Implant Stability Against Migration

The computed tomography (CT) scans of the gastrocnemius and iliotibialiscranialis muscles were used to create a migration plot (see FIG. 13 )showing the separation distances of the magnetic bead pairs over time.The magnetic bead pair that was implanted closest to one another, withan initial separation distance measured at 15.3 mm, did undergomigration, as can be seen by the 3 mm final separation of this pair inFIG. 13 . The second closest magnet pair, with an initial separationdistance measured at 16.7 mm, did not fully migrate. However, as themigration study is still ongoing at the time of this writing, this pairmay yet experience additional migration. Upon observing all othermagnetic bead pairs, it appears that the beads at longer separationdistances are all resilient long-term to migration, suggesting magneticbeads should be implanted with separation distances of at least 20 mm ormore.

Discussion: Magnetomicrometry

The results presented above suggest that magnetomicrometry is viable forhuman prosthetic and exoskeletal control. The results presented abovedemonstrate sub-millimeter accuracy, and with improved sensorpositioning or calibration, it appears that the accuracy resolution canbe easily brought to well below 100 micrometers. Even in the absence ofimproved sensor positioning, the precision of the measurements issufficient to faithfully track normalized muscle lengths and velocitiesin real time. In other words, when the geometry of an applicationpermits, it is expected that comparable to superior accuracy can beachieved from magnetomicrometry as would be found in fluoromicrometry,but with the additional benefits of being portable, low-cost, andreal-time.

The frequency sweep data that was collected was performed via passivecycling of the muscle via movement of the ankle joint, but largerexcursions are expected when the muscle is actively cycled.

The data shown in FIG. 11 , from the static measurement interferencetest, demonstrates that the process of fluoromicrometry does notinterfere with the magnetomicrometry measurements. However, theseresults do not provide insight into whether the magnetomicrometryinterferes with the fluoromicrometry.

In the single-board static magnetomicrometry experiments (see FIG. 12 ),the spread in the measurement (i.e., the standard deviation) can beexplained simply by the normally-distributed Johnson and flicker noise,but the mean error cannot be explained in this same way, so the meanoffset deserves additional investigation to determine the cause.

Discussion: Long-Term Implant Stability Against Migration

Long-term implant stability depends on the properties of muscle tissue,the size and coating of the magnets, and the forces that the magnetsexert on one another. The interaction between muscle tissue propertiesand the size and coating of the magnets on long-term stability is bestempirically investigated, but the forces that magnets exert on oneanother is given by a simple mathematical model. Using equation (37) of[24] and assuming that the magnets in each pair of magnetic beads arealigned with one another, the magnitude of the force between two magnetsin a magnet pair is given as:

$\begin{matrix}\begin{matrix}{F = {{\frac{3\mu_{0}m_{1}m_{2}}{4\pi d^{4}}\left( {{\hat{d}\left( {{\hat{m}}_{1} \cdot {\hat{m}}_{2}} \right)} + {{\hat{m}}_{1}\left( {\hat{d} \cdot {\hat{m}}_{2}} \right)} + {{\hat{m}}_{2}\left( {\hat{d} \cdot {\hat{m}}_{1}} \right)} -} \right.}}} \\{\left. {}{5{\hat{d}\left( {\hat{d} \cdot {\hat{m}}_{1}} \right)}\left( {\hat{d} \cdot {\hat{m}}_{2}} \right)} \right)} \\{{= {\frac{3\mu_{0}}{2\pi}\frac{m_{1}m_{2}}{d^{4}}}},}\end{matrix} & (1.1)\end{matrix}$

where d is the separation distance between the magnets. Thus, with theabove assumption on magnet alignment, Eqn. 1.1 can be used to determinean estimate of the force between the magnets. In Table 1.1, the magnetstrengths as measured before implantation are compiled along with thestandard deviation of these strength measurements and the minimum andmaximum measured magnet separation distances after implantation.

TABLE 1.1 Implanted Magnetic Bead Pair Strengths and SeparationsMagnetic Dipole Strength (mA · m²) Magnet 1 Magnet 2 Separation (mm)Turkey ID Side Muscle Strength Std Dev Strength Std Dev Max Min A LGastroc 10.3 2.0 10.4 0.6 15.3 3.6 ITC 9.0 0.5 9.0 0.5 21.2 19.4 RGastroc 9.1 0.5 9.2 0.4 34.6 34.0 ITC 9.2 1.9 9.3 0.5 44.2 41.5 B LGastroc 9.7 0.7 9.7 0.6 39.2 37.1 ITC 9.6 0.8 9.6 0.5 30.4 27.7 RGastroc 9.4 0.4 9.4 0.4 51.6 50.5 ITC 9.4 0.4 9.6 0.3 27.9 26.7 C LGastroc 10.2 0.5 10.2 0.5 66.5 65.6 ITC 8.1 1.5 10.1 1.4 21.5 21.0 RGastroc 9.9 0.4 9.9 0.5 69.3 68.4 ITC 10.1 0.4 9.9 0.4 16.7 14.3 D LGastroc 9.9 0.4 9.9 0.4 62.3 60.7 ITC 9.8 1.9 10.0 0.4 22.1 21.6 RGastroc 10.0 0.5 10.0 0.4 69.7 68.9 ITC 9.8 1.4 10.1 1.0 23.8 23.6

It may be tempting to use Eqn. 1.1 with Table 1.1 to calculate a maximumallowable force between magnets. For instance, using the first row ofthe Table 1.1, corresponding to the magnetic bead pair that migrated, itis tempting to consider using the strength of the migrated magnet pairminus the measurement standard deviation for each magnet, along with themaximum measured separation distance of the migrated magnet pair, todetermine a somewhat conservative estimate of maximum allowable force.This approach, however, would be flawed. Assuming a neodymium iron boronvolumetric mass density of 7.5 g/cm³, a 3 mm diameter sintered neodymiumiron boron magnet experiences a force due to gravity of approximately1.040 mN. In contrast, the above calculation gives only 0.887 mN,significantly less than the force due to gravity on the magnet.Revisiting Eqn. 1.1, the force increases with the inverse fourth powerof the distance between the magnetic beads, so, importantly, themagnetic beads experience a significantly greater force at smallerdistances. For instance, now using just the magnetic dipole strengthsfrom the first row of Table 1.1 with Eqn. 1.1, the force versus magneticbead separation distance is shown in FIG. 14 .

Noting this steep increase in force as the magnetic bead separationdistance is reduced, it is essential to note that the minimum magneticbead separation distance should take into account the minimum distancethat the magnetic beads will come within one another under maximumcontraction. As seen in the magnetomicrometry studies, the magnets comeseveral millimeters closer to one another during passive cycling. It isexpected that the magnetic beads come even closer to one another duringactive flexion of the muscles.

Finally, migration is also dependent upon both the magnets' strength andthe magnets' geometry, so one can take into account tolerances in bothof these parameters when creating specifications for safe magnetseparation distances in tissue.

BIBLIOGRAPHY

1. Ariel L Camp, Henry C Astley, Angela M Horner, Thomas J Roberts, andElizabeth L Brainerd. Fluoromicrometry: a method for measuring musclelength dynamics with biplanar videofluoroscopy. Journal of ExperimentalZoology Part A: Ecological Genetics and Physiology, 325(7):399-408,2016.

2. Jared W Coburn, T J Housh, J T Cramer, J P Weir, J M Miller, T WBeck, M H Malek, and G O Johnson. Mechanomyographic time and frequencydomain responses of the vastus medialis muscle during submaximal tomaximal isometric and isokinetic muscle actions. Electromyography andclinical neurophysiology, 44(4):247-255, 2004.

3. John V Frangioni, Tao S Kwan-Gett, Lynn E Dobrunz, and Thomas AMcMahon. The mechanism of low-frequency sound production in muscle.Biophysical journal, 51(5):775-783, 1987.

4. Wolfgang Kabsch. A solution for the best rotation to relate two setsof vectors. Acta Crystallographica Section A: Crystal Physics,Diffraction, Theoretical and General Crystallography, 32(5):922-923,1976.

5. Wolfgang Kabsch. A discussion of the solution for the best rotationto relate two sets of vectors. Acta Crystallographica Section A: CrystalPhysics, Diffraction, Theoretical and General Crystallography,34(5):827-828, 1978.

6. K&J Magnetics. Neodymium magnet physical properties: Magnet summarytable.http://web.archive.org/web/20200331161055/https://www.kjmagnetics.com/specs. asp. Accessed: 2020-03-31.

7. SM Magnetics. Magnetic materials & tables: Neodymium.https://smmagnetics. com/pages/magnetic-materials-tables. Accessed: 2020Mar. 31, Archived athttp://web.archive.org/web/20200331160201/https://cdn.shopify.com/s/files/1/1359/3085/files/Material_tables_on_website_-Neodymium_-_2.pdf?4920142403209868357.

8. Thomas A McMahon. “Chapter 1: Fundamental Muscle Mechanics”. Muscles,reflexes, and locomotion. Princeton University Press, 1984.

9. T Nara and W Ito. Moore—penrose generalized inverse of the gradienttensor in euler's equation for locating a magnetic dipole. Journal ofApplied Physics, 115(17):17E504, 2014.

10. Takaaki Nara, Satoshi Suzuki, and Shigeru Ando. A closed-formformula for magnetic dipole localization by measurement of its magneticfield and spatial gradients. IEEE Transactions on Magnetics,42(10):3291-3293, 2006.

11. Talat Ozyagcilar. Calibrating an eCompass in the presence of hardand soft-iron interference. Freescale Semiconductor Ltd, pages 1-17,2012.

12. None

13. Timo Pylvanainen. Automatic and adaptive calibration of 3d fieldsensors. ^({umlaut over ( )}) Applied Mathematical Modelling,32(4):575-587, 2008.

14. Thomas J Roberts, Richard L Marsh, Peter G Weyand, and C RichardTaylor. Muscular force in running turkeys: the economy of minimizingwork. Science, 275(5303):1113-1115, 1997.

15. Craig Schroeder. Practical course on computing derivatives in code.In ACM SIGGRAPH 2019 Courses, pages 1-22. 2019.

16. S Tarantino, F Clemente, D Barone, M Controzzi, and C J S RCipriani. The myokinetic control interface: tracking implanted magnetsas a means for prosthetic control. Scientific reports, 7(1):1-11, 2017.

17. Cameron R Taylor, Haley G Abramson, and Hugh M Herr. Low-latencytracking of multiple permanent magnets. IEEE Sensors Journal,19(23):11458-11468, 2019.

18. Arnold Magnetic Technologies. Neodymium iron boron magnet catalog:Sintered neodymium-iron-boron magnets; other properties.https://web.archive.org/web/20200513225014/https://www.arnoldmagnetics.com/wp-content/uploads/2019/06/Arnold-Neo-Catalog.pdf.Accessed: 2020 May 13.

19. Arnold Magnetic Technologies. Neodymium iron boron magnets:Available neo grades.http://web.archive.org/web/20200331161416/https://www.arnoldmagnetics.com/products/neodymium-iron-boron-magnets/. Accessed: 2020 Mar. 31.

20. D A Turner, I J Anderson, J C Mason, and M G Cox. An algorithm forfitting an ellipsoid to data. National Physical Laboratory, UK, 1999.

21. The UNIX and Linux Forums. Linux and unix man pages: Bsd libraryfunctions manual; sincos(3).https://www.unix.com/man-page/osx/3/_SINCOS/. Accessed: 2020 May 18.

22. R F Wiegert. Magnetic star technology for real-time localization andclassification of unexploded ordnance and buried mines. In Detection andSensing of Mines, Explosive Objects, and Obscured Targets XIV, volume7303, page 73031U. International Society for Optics and Photonics, 2009.

23. W Wynn, C Frahm, P Carroll, R Clark, J Wellhoner, and M Wynn.Advanced superconducting gradiometer/magnetometer arrays and a novelsignal processing technique. IEEE Transactions on Magnetics,11(2):701-707, 1975.

24. Kar W Yung, Peter B Landecker, and Daniel D Villani. An analyticsolution for the force between two magnetic dipoles. Physical Separationin Science and Engineering, 9(1):39-52, 1998.

25. Kathryn Ziegler-Graham, Ellen J MacKenzie, Patti L Ephraim, Thomas GTravison, and Ron Brookmeyer. Estimating the prevalence of limb loss inthe united states: 2005 to 2050. Archives of physical medicine andrehabilitation, 89(3):422-429, 2008.

26. Martin, Jack A., et al. “Gauging force by tapping tendons.” Naturecommunications 9.1 (2018): 1-9.

The teachings of all patents, published applications and referencescited herein are incorporated by reference in their entirety.

While example embodiments have been particularly shown and described, itwill be understood by those skilled in the art that various changes inform and details may be made therein without departing from the scope ofthe embodiments encompassed by the appended claims.

What is claimed is:
 1. A method of detecting muscle activation,comprising: with a magnetic field sensor, detecting lateral vibration ofa target implanted at a muscle or a tendon, the target comprising amagnetic material; and estimating a level of muscle activation based onthe detected lateral vibration.
 2. The method of claim 1, furthercomprising estimating a muscle force based on the estimated level ofmuscle activation and a length and a velocity of the muscle.
 3. Themethod of claim 1 or 2, wherein detecting lateral vibration includesdetecting movements of the target in a range of about 10 μm to about 20μm.
 4. The method of any one of claims 1-3, wherein detecting lateralvibration includes detecting vibrational movement of the target atfrequencies of greater than about 10 Hz.
 5. The method of any one ofclaims 1-4, further comprising, with an accelerometer, detectingvibrational movement of the magnetic field sensor relative to thetarget.
 6. The method of any one of claims 1-5, wherein detectinglateral vibration of a target includes detecting lateral vibrations oftwo or more targets disposed along an axis.
 7. A method of estimating aphysiological parameter of a muscle or tendon, comprising: with amagnetic field sensor, detecting vibration of a target implanted at amuscle or tendon, the target comprising a magnetic material; andestimating at least one of a muscle force and tendon force based on thedetected vibration.
 8. The method of claim 7, further comprising:applying a perturbance to the target or to a muscle-tendon unitcomprising the muscle and the tendon; measuring a timing of thevibration of the target relative to a timing of the perturbance orrelative to a timing of a vibration of one or more additional targetsimplanted at the muscle-tendon unit; and estimating a speed of a shearwave or a compression wave in the muscle-tendon unit or surroundingtissue of the muscle-tendon unit based on the measured timing.
 9. Themethod of claim 8, wherein applying a perturbance includes poking orvibrating the muscle-tendon unit or the surrounding tissue.
 10. Themethod of claim 9, wherein the poking or vibrating includes actuating amagnetic bead affixed at the muscle tendon unit or the surroundingtissue by an applied magnetic field.
 11. The method of claim 10, whereinthe applied magnetic field is supplied by an electromagnet.
 12. Themethod of claim 11, wherein the electromagnet is external to the body.13. The method of claim 8, further comprising determining aphysiological property of the muscle-tendon unit based on the estimatedspeed of the shear wave or the compression wave.
 14. The method of claim13, wherein the physiological property is stiffness of the muscle, thetendon, or the surrounding tissue.
 15. A method of monitoringbiomechanical motion, comprising: disposing at least one target on asubject at a location associated with a joint of the subject; with amagnetic field sensor array, detecting a change in state of the at leastone target relative to the magnetic field sensor array, another targetdisposed on the subject, or a combination thereof; and determining astate of the joint based on the detected change in state of the at leastone target.
 16. A method for determining one or more of three sensorposition parameters and three sensor orientation parameters for each ofthe sensors in a sensor array, comprising: placing at least one targetin at least one known location relative to a sensor array, whereby asignal from the at least one target at the sensors is detected, andrecording at least one measurement of the signal at each of the sensorsfor each placement of the one or more targets; estimating one or moreparameters from the group consisting of x-position, y-position,z-position, yaw, pitch, and roll, of each of the sensors; estimating anyunknown state parameters of the at least one target; providing aconstant value for a magnetic dipole weight state parameter of the atleast one target for each of the measurements; calculating predictedvalues of the signal at each of the sensors for each of the measurementsgiven the one or more estimated sensor parameters, the estimated targetstate parameters, and the provided constant magnetic dipole weight stateparameter; computing a prediction error in the predicted values of thesignal with reference to the values of the signals detected at thesensors; calculating a prediction error Jacobian matrix by analyticallycomputing elements of the prediction error Jacobian matrix with respectto the estimated parameters of the sensors for each measurement; anddetermining from the prediction error and the prediction error Jacobianmatrix a state of the parameters of the sensors.
 17. A method ofcalibrating a magnetic field sensor array, comprising:three-dimensionally rotating a magnetic field sensor array in a uniformmagnetic field, the magnetic field sensor array comprising a pluralityof sensors; recording data from each of the plurality of sensors;calculating a non-rotating transformation of the recorded data using anellipsoid fit of the data for each of the plurality of sensors;calculating a scaling factor of each of the plurality of sensors basedon relative dimensions of the ellipsoid fit; applying the calculatednon-rotating transformation and calculated scaling factor to obtain atransformed, scaled dataset for each of the plurality of sensors; androtating the obtained transformed, scaled datasets to align, therebydetermining a relative orientation of each of the plurality of sensorsto calibrate the magnetic field sensor array.
 18. An implantable targetfor magnetic tracking, comprising: a base structure comprising amagnetic or magnetizable material; and a shell structure disposed aboutthe base structure and comprising layers of nickel, copper, gold, andparylene C.
 19. The implantable target of claim 18, wherein the layersof the shell structure are arranged from the base structure in order ofnickel, copper, nickel, gold, and parylene C.
 20. The implantable targetof claim 18 or 19, wherein the gold layer of the shell structure has athickness of at least about 5 μm.
 21. The implantable target of any oneof claims 18-20, wherein the parylene C layer of the shell structure hasa thickness of at least about 25 μm.
 22. The implantable target of anyone of claims 18-21, wherein the magnetic or magnetizable materialcomprises neodymium iron boron and dysprosium.
 23. An insertion devicefor an implantable target, comprising: a cannula; a cartridge configuredto position an implantable target at the cannula; and a pushrodreceivable in the cannula and configured to push the implantable targetfrom the cartridge through the cannula for delivery to a tissue, thecannula and the pushrod comprising a nonmagnetic material.
 24. Theinsertion device of claim 23, wherein a distal end of the pushrod is ofa complementary geometry to the implantable target to prevent damage toa shell structure of the target during delivery.
 25. The insertiondevice of claim 23, wherein a distal surface of the pushrod comprises aspring structure to prevent damage to a shell structure of the targetduring delivery.
 26. The insertion device of any one of claims 23-25,further comprising a mount coupling the cartridge to the cannula,wherein the cartridge is configured to retain a plurality of targets andis rotatable about the mount to position each of the plurality oftargets at the cannula.
 27. An insertion system, comprising: a pluralityof the implantable targets of claim 18; and the insertion device ofclaim 23, the plurality of implantable targets being deliverable by theinsertion device.
 28. A wearable shielding assembly, comprising: anarray of sensors configured to detect a state change of at least onemagnetic target implanted at a tissue; a wearable receptacle withinwhich the array of sensors is disposed; and ageometrically-reconfigurable material disposed about or integral withthe wearable receptacle and configured to provide magnetic shielding tothe array of sensors.
 29. The assembly of claim 28, wherein thegeometrically-reconfigurable material is flexible.
 30. The assembly ofclaim 28, wherein the geometrically-reconfigurable material comprises aplurality of ferromagnetic tiles.
 31. The assembly of claim 30, whereinthe ferromagnetic tiles are releasably connectable to one another. 32.The assembly of claim 28, wherein the geometrically-reconfigurablematerial comprises a ferromagnetic mesh.
 33. The assembly of claim 28,wherein the geometrically-reconfigurable material is rigidly affixedafter application to the wearable receptacle.
 34. A method of trackingone or more permanent magnets, comprising: detecting a signal from eachof the one or more permanent magnets; calculating ananalytically-derived Hessian matrix with respect to the detectedsignals; and determining a state of each of the one or more permanentmagnets based on the calculated Hessian matrix.
 35. The method of claim34, further comprising: calculating a predicted magnetic field value foreach of the one or more permanent magnets; and calculating a Jacobianmatrix with respect to the detected signals, wherein determining thestate of each of the one or more permanent magnets is further based onthe calculated predicted magnetic field values and the calculatedJacobian matrix, and wherein calculation of the Hessian matrix, theJacobian matrix, and the predicted magnetic field values is performedsimultaneously using common subexpression elimination.
 36. The method ofclaim 34, wherein the calculation of the Hessian matrix is parallelized.37. A method of assembling a non-planar sensing array, comprising:fabricating a plurality of sensors on a flexible circuit board; andaffixing the flexible circuit board to a non-planar and rigid substrate.38. The method of claim 37, wherein the substrate is a prostheticsocket.
 39. The method of claim 37 or 38, wherein the sensors aremagnetic field sensors.
 40. The method of any one of claims 37-39,wherein the flexible circuit board is a long strip.
 41. A system,comprising: at least two sensor arrays, each sensor array configured todetect a state of at least one magnetic target at a tissue; and at leastone position sensor associated with at least one of the at least twosensor arrays and configured to detect a position and orientation of theassociated sensor array relative to the other of the at least two sensorarrays.
 42. The system of claim 41, wherein the at least one positionsensor is an inertial measurement unit disposed at the associated sensorarray.
 43. The system of claim 41 or 42, further comprising a controllerconfigured to: determine a distance and difference in orientationbetween each of the at least two sensor arrays and the at least onemagnetic target; and associate the at least one magnetic target with oneof the at least two sensor arrays based upon the determined distancesand orientations.
 44. The system of any one of claims 41-43, wherein theat least one magnetic target is implanted at the tissue.
 45. A systemfor detecting muscle activation, comprising: a magnetic field sensorconfigured to detect lateral vibration of at least one target implantedat a muscle or a tendon, the at least one target comprising a magneticmaterial; and a controller configured to estimate a level of muscleactivation based on the detected lateral vibration.
 46. The system ofclaim 45, wherein the controller is further configured to estimate amuscle force based on the estimated level of muscle activation and alength and a velocity of the muscle.
 47. The system of claim 45 or claim46, further comprising an accelerometer at the magnetic field sensor andconfigured to detect vibrational movement of the magnetic field sensorrelative to the target.
 48. The system of any one of claims 45-47,wherein the at least one target comprises two or more targets disposedalong an axis.
 49. A system for estimating a physiological parameter ofa muscle or tendon, comprising: a magnetic field sensor configured todetect vibration of at least one target implanted at a muscle or atendon, the at least one target comprising a magnetic material; and acontroller configured to estimate at least one of a muscle force, atendon force, and a muscle-tendon unit force based on the detectedvibration.
 50. The system of claim 49, wherein the controller is furtherconfigured to estimate a speed of a shear wave or a compression wave inthe muscle, the tendon, a muscle-tendon unit comprising the muscle andthe tendon, or surrounding tissue of the muscle-tendon unit based on ameasured timing of a vibration of the at least one target relative to atiming of a perturbance applied to the target or relative to one or moreadditional targets implanted at the muscle or tendon.